Dear Prof. Ben Cox and Prof. Bradley Treeby, Dear colleagues,
I am doing 2D image reconstruction with point sensor ring array using 'universal back-projection algorithm (UBP). As you know, the UBP in Prof. Lihong Wang's article 'Universal back-projection algorithm for photoacoustic computed tomography' was obtained in 3D, so I did a simplification of his reconstruction equation from 3D to 2D (To be honest, I am not 100% sure if I did it correctly. But here on the forum, I cannot post equations and pics).
I used k-wave tool box to generate the forward date, from which I did the image reconstruction. However, my UBP reconstruction has some wired fluctuations at the image transition parts compared to time reversal image reconstruction, which I really don't know the reason. Increasing the number of sensors and decreasing the grid size do not alleviate the above mention reconstruction fluctuations much. I am running out thoughts now!
Below, I post my matlab code. Could anyone please have a look at my code to see where I may not do correctly?
If anyone has done 2D universal back-projection image reconstruction successfully, could you please kindly share me how you did it？
In case you need further input, my email: isabelle.thomas19950321@gmail.com.
Thank you all!
Best wishes,
Isabelle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MATLAB CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars
close all
% load the initial pressure distribution as shepp-Logan phantom
p0_magnitude = 1;
grid_size_SL = 100;
p = p0_magnitude.*phantom('Modified Shepp-Logan',grid_size_SL);
tot_grid_size = grid_size_SL + 100; %% the extra 100 grid points are used for PML and placing transducers
p0 = zeros(tot_grid_size);
p0(51:tot_grid_size-50,51:tot_grid_size-50) = p;
% assign the grid size and create the computational grid
PML_size = 20; % size of the PML in grid points
Nx = tot_grid_size; % number of grid points in the x direction
Ny = tot_grid_size; % number of grid points in the y direction
x= 10e-3; % total grid size [m]
y = 10e-3; % total grid size [m]
dx = x / Nx; % grid point spacing in the x direction [m]
dy = y / Ny; % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% smooth the initial pressure distribution and restore the magnitude
p0 = smooth(kgrid, p0, true);
% assign to the source structure
source.p0 = p0;
% define the properties of the propagation medium
sound_speed = 1500;
medium.sound_speed = sound_speed; % [m/s]
% define a centered Cartesian circular sensor
sensor_radius = 4.5e-3; % [m]
sensor_angle = 2*pi; % [rad]
sensor_pos = [0, 0]; % [m]
num_sensor_points = 100;
cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle);
% assign to sensor structure
sensor.mask = cart_sensor_mask;
% create the time array
kgrid.makeTime(medium.sound_speed);
% set the input options
input_args = {'Smooth', false, 'PMLInside', true, 'PlotPML', true};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
%% back projection
t_array = kgrid.t_array;
dt = kgrid.dt;
fs = 1/dt;
np = length(t_array);
NFFT = 2^nextpow2(np);
f_vec = fs*linspace(0,1,NFFT);
w_vec = 2*pi*f_vec;
k_vec = ifftshift(2*pi*f_vec/(sound_speed));
weighting_omega = 0;
p0_recon = zeros(Nx, Ny);
for ii = 1:num_sensor_points
x_ii = cart_sensor_mask(1,ii);
y_ii = cart_sensor_mask(2,ii);
distance_x_square = (kgrid.x - x_ii).^2;
distance_y_square = (kgrid.y - y_ii).^2;
distance_xy = sqrt(distance_x_square + distance_y_square);
distance_xy_index = round(distance_xy/(sound_speed*dt));
distance_xy_index(~distance_xy_index)=1;
p_ii = sensor_data(ii,:);
p_ii_derivative0 = ifft((1i*k_vec).*fft(p_ii,NFFT));
p_ii_derivative = real((p_ii_derivative0(1:np)));
bp_ii = (p_ii - sound_speed.*t_array.*p_ii_derivative);
p0_ii = bp_ii(distance_xy_index);
%%% weighting angle
weighting_arc_ii = (2*pi*sensor_radius)/num_sensor_points;
weighting_phy_ii = 2*pi/num_sensor_points;
weighting_omega_ii = weighting_arc_ii*(1./((distance_x_square + distance_y_square).*sensor_radius).*(-x_ii.*kgrid.x - y_ii.*kgrid.y + sensor_radius.^2));
weighting_omega_ii(isnan(weighting_omega_ii)) = weighting_phy_ii;
weighting_omega_ii(isinf(weighting_omega_ii)) = weighting_phy_ii;
p0_ii = weighting_omega_ii.*p0_ii;
weighting_omega = weighting_omega + weighting_omega_ii;
p0_recon = p0_recon + p0_ii;
end
p0_recon0 = p0_recon./(2*pi);
figure(1) %% reconstructed image
imagesc(kgrid.y_vec*1000,kgrid.x_vec*1000,p0_recon0,[-1,1]) %./max(p0_recon(:))
axis image
colorbar
ylabel('x(mm)');
xlabel('y(mm)');
slice_pos = 4.5e-3; %% comparison with original image amplitude along a horizontal line
figure(2)
plot(kgrid.y_vec,p0(round(slice_pos/kgrid.dx),:))
hold on
plot(kgrid.y_vec,p0_recon0(round(slice_pos/kgrid.dx),:))
hold off
axis tight
xlabel('y(mm)');
ylabel('Reconstructed amplitude');