Hi,

I reconstructed initial pressure using 3D Time Reversal Reconstruction For A Planar Sensor in heterogeneous medium. But i can not get good reconstructed images.The codes I use are as following:

%% create the computational grid

PML_size = 10; % size of the PML in grid points

% Nx = 64 * scale - 2 * PML_size; % number of grid points in the x direction

% Ny = 64 * scale - 2 * PML_size; % number of grid points in the y direction

% Nz = 32 * scale - 2 * PML_size; % number of grid points in the z direction

Nx = 190 - 2 * PML_size; % number of grid points in the x direction

Ny = 270 - 2 * PML_size; % number of grid points in the y direction

Nz = 126 - 2 * PML_size; % number of grid points in the z direction

dx = 1.25e-3/2 ; % grid point spacing in the x direction [m]

dy = 1.25e-3/2; % grid point spacing in the y direction [m]

dz = 1.25e-3/2; % grid point spacing in the z direction [m]

kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

%% load the speed of sound and the tissue density and the initial pressure distribution from an image and scale

soundspeed=importdata('C:\Users\WANG MENG XIAO\Desktop\project\XACT\code_and_results\work3_06152018\a_segmentation\mat\soundspeed_160_240_96.mat');

tissdst=importdata('C:\Users\WANG MENG XIAO\Desktop\project\XACT\code_and_results\work3_06152018\a_segmentation\mat\tissdst_160_240_96.mat');

p0=importdata('C:\Users\WANG MENG XIAO\Desktop\project\XACT\code_and_results\work3_06152018\b_initial_pressure\mat\initial_dose_160_240_53.mat');

%% get the size of soundspeed

[r,c,h]=size(soundspeed);

%% define the speed of sound(m/s)

for k=1:h

for i=1:r

for j=1:c

if soundspeed(i,j,k)<=343

soundspeed(i,j,k)=1500;

end

end

end

end

%% define the tissue density(kg/m3)

for k=1:h

for i=1:r

for j=1:c

if tissdst(i,j,k)<=1.2

tissdst(i,j,k)=1000;

end

end

end

end

p0_binary=zeros(Nx,Ny,Nz);

p0_binary(6:165,6:245,6:58)=p0;

% p0_binary=zeros(Nx,Ny,Nz);

% p0_binary(1:32,1:32,5)=im;

% smooth the initial pressure distribution and restore the magnitude

% p0 = smooth(kgrid, p0_binary, true);

p0=p0_binary;

%% assign to the source structure

source.p0 = p0;

%% define the properties of the propagation medium

sound_speed=ones(Nx,Ny,Nz)*1500;

sound_speed(6:165,6:245,6:101)=soundspeed;

medium.sound_speed = sound_speed; % [m/s]

density=ones(Nx,Ny,Nz)*1000;

density(6:165,6:245,6:101)=tissdst;

medium.density=density; % [kg/m3]

% 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(6:2:165,6:2:245,102) = 1;

% sensor.mask(1:32, 1:32,1) = 1;

%% define frequence and bandwidth of sensor

center_freq =1e6; %!!!!!!定义探测器的中心频率

bandwidth = 100; %带宽100%

sensor.frequency_response = [center_freq, bandwidth]; %定义传感器完整参数，包括中心频率和带宽

%% create the time array

kgrid.makeTime(medium.sound_speed);

%% set the input arguements

input_args = {'PMLSize', PML_size, 'PMLInside', false, ...

'PlotPML', false, 'Smooth', false, 'DataCast', 'single'};

%% run the simulation

sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

%% reset the initial pressure

source.p0 = 0;

%% assign the time reversal data

sensor.time_reversal_boundary_data = sensor_data;

%% run the time-reversal reconstruction

p0_recon = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});

%% add first order compensation for only recording over a half plane

p0_recon = 2 * p0_recon;

%% apply a positivity condition

p0_recon(p0_recon < 0) = 0;

Thank you very much.

Sincerely,

Mengxiao