3D reconstruction using a plane transducer

Hi dear bencox,
Here is my code:
%% define literals
PML_size = 10; % size of the perfectly matched layer in grid points

%% create the computational grid
Nx = 120 - 2 * PML_size; % number of grid points in the x direction
Ny = 120 - 2 * PML_size; % number of grid points in the y direction
Nz = 120 - 2 * PML_size; % number of grid points in the z direction
dx = 1.25e-3 ; % grid point spacing in the x direction [m]
dy = 1.25e-3; % grid point spacing in the y direction [m]
dz = 1.25e-3; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

%% load the initial pressure distribution from an image and scale
p0=importdata('C:\Users\WANG MENG XIAO\Desktop\TR_reconstruction\codes\initial_pressure.mat');
p0_binary=zeros(Nx,Ny,Nz);
p0_binary(39:62,35:66,26:43)=p0;

p0=p0_binary;

%% assign to the source structure
source.p0 = p0;

%% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
medium.density=1000; % [kg/m^3]
medium.alpha_coeff = 0.75; % power law absorption prefactor [dB/(MHz^y*cm)]
medium.alpha_power = 1.5; % power law absorption exponent

%% define a binary planar sensor
sensor.mask = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz);
sensor.mask(35:66,35:66,76) = 1;

%% define frequence and bandwidth of sensor
center_freq =1e6; % the center frequency of sensor(Hz)
bandwidth = 100; % bandwidth 100%
sensor.frequency_response = [center_freq, bandwidth];

%% create the time array
t_end=sqrt((Nx*dx).^2 + (Ny*dy).^2 + (Nz*dz).^2)/medium.sound_speed;
cfl = 0.3;
dt = cfl*dx/medium.sound_speed;
kgrid.t_array = 0:dt:t_end;

%% set the input arguements
input_args = {'PMLSize', PML_size, 'PMLInside', false, ...
'PlotPML', false, 'Smooth', false, 'PlotSim', true};

%% run the simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

%% =========================================================================
% RECONSTRUCT AN IMAGE USING TIME REVERSAL
% =========================================================================

% remove the initial pressure field used in the simulation
source = rmfield(source, 'p0');

% use the sensor points as sources in time reversal
source.p_mask = sensor.mask;

%% time reverse and assign the data
%source.p=zeros(sum(source.p_mask(:)),kgrid.Nt);
source.p = fliplr(sensor_data);

% enforce, rather than add, the time-reversed pressure values
source.p_mode = 'dirichlet';

% set the simulation to record the final image (at t = 0)
sensor.record = {'p','p_final'};

% run the time reversal reconstruction
p0_estimate = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

%% apply a positivity condition
% p0_estimate.p_final = p0_estimate.p_final .* (p0_estimate.p_final > 0);

% store the latest image estimate
p0_1 = p0_estimate.p_final;

The part that is close to and farther awway from the detectors have good reconstruction results. However, the middle part has bad results.
I don't know how to post images. If you dont mind, would you tell me your email address so I can sent it to you?

k-Wave

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields