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)