Dear developers,
Thank you so much for providing this wonderful simulator!
We are trying to generate RF data for photoacoustic computed tomography that uses multiple linear arrays circularly positioned around the target.
We tried the code below and performed time reversal for image reconstruction.
However, we got a wrong reconstruction that looks like the time reversing was done a little bit more.
(i.e. circular targets appearing as donut shapes)
We then tried to reconstruct the image using self-developed delay-and-sum algorithm, and we still got a similar result with time reversal.
We found an unexpected error in the RF data: the wavefronts appearing at earlier timings than expected.
Using point source targets, we calculated the expected timing of the wavefront, and we found that the wavefronts in the RF data appeared at about 88% of the calculated delay.
When we manually multiplied 0.88 at the calculated delays for delay-and-sum, the image was well reconstructed.
We think time reversal operation was working all right, but the forward process was somehow not right.
We have provided the code as the following lines.
Sorry that we couldn't include the element positions here for you to test the code.
If you need to execute the code, you can simply put a single 128-element linear array on the top of the grid.
(We also tried that but it also gave the same error.)
Thank you in advance for kindly taking your time on this issue.
Wonseok Choi.
% Environment setting
PML_size = 40; % size of the PML in grid points --> PML = perfectly matched layer. assumed to simulate open boundaries
Nz = 1600; % number of grid points in the x (row) direction
Nx = 1600; % number of grid points in the y (column) direction
dz = 1e-4; % grid point spacing in the z direction [m]
dx = 1e-4; % grid point spacing in the x direction [m]
Xsize = Nx * dx; Zsize = Nz*dz;
% Detector position
load('ElementPosition.mat')
% define a binary line sensor
sensor.mask = zeros(Nz, Nx);
posX = pos_ele(:,1);
posXgrid = round((posX-(max(posX)+min(posX))/2)/abs(dx)) + Nx/2;
posZ = pos_ele(:,2);
posZgrid = round((posZ-(max(posZ)+min(posZ))/2)/abs(dz)) + Nz/2;
for idxEle = 1:length(posX)/8
sensor.mask(posZgrid(idxEle), posXgrid(idxEle)) = 1;
end
% Temporal arguments
SoundSpeed = 1510;
SampleFreq = 50e6;
dt = 1/SampleFreq;
Nt = ceil(1.1 * sqrt(Xsize^2+Zsize^2)/SoundSpeed/dt)/2;
% Transducer arguments
transducer.pitch = 0.298e-3; % [m]
transducer.numElement = 128; % [m]
transducer.centerFreq = 5e6; % [Hz]
transducer.BW = 60; % [%]
sensor.frequency_response = [transducer.centerFreq, transducer.BW];
% create the computational grid
kgrid = kWaveGrid(Nz, dz, Nx, dx);
kgrid.setTime(Nt, dt);
% Medium Parameters
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.01;
medium.sound_speed = SoundSpeed; % [m/s]
% create initial pressure distribution using makeDisc
NumSource_x = 3;
NumSource_z = 3;
Disc.magnitude = 1; % [Pa]
Disc.radius = 2; % scaled by grid point spacing.
numTarget = 9;
disc = 0;
for idx_x = 1:NumSource_x
for idx_z = 1:NumSource_z
Disc.z_pos = 384*(idx_z-2)/(NumSource_z+1)+Nz/2; % [grid points]
Disc.x_pos = 384*(idx_x-2)/(NumSource_x+1)+Nx/2; % [grid points]
temp = Disc.magnitude * makeDisc(Nz, Nx, Disc.z_pos, Disc.x_pos, Disc.radius);
disc = disc + temp;
end
end
% smooth the initial pressure distribution and restore the magnitude
source.p0 = smooth(kgrid, disc, true);
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false};
sensor_data_raw = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% assign the time reversal data
source2 = source;
source2.p0=0;
sensor2 = sensor;
sensor2.time_reversal_boundary_data = sensor_data_raw;
kgrid2 = kgrid;
medium2 = medium;
medium2.alpha_coeff = 0;
p0_tr = kspaceFirstOrder2D(kgrid2, medium2, source2, sensor2, input_args{:});