Hi,

I am trying to generate a basic TR reconstruction (3d) for planar geometry. Although, my problem involves reconstructing metal-tissue interface, I am at first trying to get reconstructions for homogenous domain settings.

I am getting a very poor reconstruction for a data-rich problem. I do expect limited view artifacts, but this is much worse than that. Please let me know if I am doing something wrong here? How can I improve my reconstructions?

Also, while carrying out reconstruction for heterogeneous case (sound speed and density can be upto 4-6 times of that of water) can I expect a good reconstruction?? Please give me insights on this aspect as well.

Sincerely,

Kevin

The problem setting, true PA source, and reconstruction are uploaded at:

https://imgur.com/a/XGRelbN

%% ************ Code TR - homogeneous SOS and density ********************

% =========================================================================

% SIMULATION

% =========================================================================

% % assign the grid size and create the computational grid

domain_size=2e-2;

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

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

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

dz =0.1e-3; % grid point spacing in the y direction [m]

Nx = round(domain_size/dx); % number of grid points in the x direction

Ny = round(domain_size/dx/5) ; %number of grid points in the y direction

Nz = round(domain_size/dx); % number of grid points in the z direction

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

p0 = zeros(Nx,Ny,Nz);

medium.sound_speed = 1500*ones(Nx,Ny,Nz); %% background

medium.density = 1000*ones(Nx,Ny,Nz); %% background

metal_sze = 2e-3; % side length of the metal cube

sz_vox = round(metal_sze/dx);

tissue_sze = 1e-3;

tissue_vox = tissue_sze/dx;

for i = 1:Nx

for j=1:Ny

for k=1:Nz

if abs(Nx/2-i)<= sz_vox/2 && abs(Ny/2-j)<=sz_vox/2 && abs(Nz/2-k)<=sz_vox/2

p0(i,j,k)=2;

% medium.sound_speed(i,j,k)=1500; %6000;

% medium.density(i,j,k)= 1000; %3000;

elseif abs(Nx/2+sz_vox/2-i)<=tissue_vox && abs(Ny/2-j)<=sz_vox/2 && abs(Nz/2-k)<=sz_vox/2

p0(i,j,k)=1;

% medium.sound_speed(i,j,k)=1500; %1600;

% medium.density(i,j,k)= 1000; %1500;

end

end

end

end

% smooth the initial pressure distribution and restore the magnitude

source.p0 = smooth(kgrid, p0, true);

% medium.sound_speed = smooth(kgrid, medium.sound_speed, true);

% medium.density = smooth(kgrid, medium.density, true);

p_slice = p0(:,:,round(Nz/2));

figure, imagesc(p_slice)

% Detection setting and measurement generation :----

% define a binary planar sensor

sensor.mask = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz);

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

%% create the time array

% [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

dt = 1/8e6;

kgrid.t_array= 0:dt:100*dt;

% set the input arguements

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

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

% run the simulation

tic

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

toc

% 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;

p_rec_slice = p0_recon(:,:,round(Nz/2));

figure, imagesc(p_rec_slice)