Hi,
I noticed that the time-reversal example code 'example_pr_2D_TR_circular_sensor.m' yields quantitatively accurate results for high number of detectors placed around the full circle.
So, I tried to see if the same can be obtained if we temporally flip the signals and send them from the transducers as a temporally varying source.
However, in this though the shape of the blood vessel is well reconstructed, the reconstructed values are about one-third of what I expected. I was wondering if someone could help me understand the reason or let me know what I am doing wrong (both codes are appended below).
The figure showing both reconstruction:
https://ibb.co/CvYjFgK
Thanks,
Kevin
[Please ignore the inverse crime -- I am running the forward and backward simulations on the same grid.]
******* Code 1 ******modified TR example****
clear
% load the initial pressure distribution from an image and scale
p0 = loadImage('EXAMPLE_source_two.bmp');
PML_size = 20; % size of the PML in grid points
Nx = 256 - 2 * PML_size; % number of grid points in the x direction
Ny = 256 - 2 * PML_size; % number of grid points in the y direction
x = 11e-2; % total grid size [m]
y = 11e-2; % 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);
p0 = resize(p0, [Nx, Ny]);
p0 = smooth(p0, true);
source.p0 = p0;
medium.sound_speed = 1500; % [m/s]
sensor_radius = 5e-2; % [m]
sensor_angle = 2*pi; % [rad]
sensor_pos = [0, 0]; % [m]
num_sensor_points = 400;
cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle);
sensor.mask = cart_sensor_mask;
kgrid.dt = 5e-8;
kgrid.t_array = [0:1500]*kgrid.dt;
input_args = {'Smooth', false, 'PMLInside', false, 'PlotPML', false};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
%****TR based Recn***
kgrid_recon = kWaveGrid(Nx, dx, Ny, dy);
kgrid_recon.setTime(kgrid.Nt, kgrid.dt);
source.p0 = 0;
sensor.time_reversal_boundary_data = sensor_data;
p0_recon = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:});
figure, imagesc(p0_recon),colormap hot,colorbar
******* Code 2 ******reconstruction by simulating detectors as a time varying source****
close all
clear all
clear
clc
% =========================================================================
% SIMULATION
% =========================================================================
% load the initial pressure distribution from an image and scale
p0_magnitude = 1; % bar
p0 = p0_magnitude * loadImage('EXAMPLE_source_two.bmp');
PML_size = 20; % size of the PML in grid points
Nx = 2*256 - 2 * PML_size; % number of grid points in the x direction
Ny = 2*256 - 2 * PML_size; % number of grid points in the y direction
x = 11e-2; % total grid size [m]
y = 11e-2; % 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);
p0 = resize(p0, [Nx, Ny]);
p0 = smooth(p0, true);
source.p0 = p0;
medium.sound_speed = 1500; % [m/s]
sensor_radius = 5e-2; % [m]
sensor_angle = 2*pi; % [rad]
sensor_pos = [0, 0]; % [m]
num_sensor_points = 400;
cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle);
[grid_data, order_index, reorder_index] = cart2grid(kgrid, cart_sensor_mask);
sensor.mask = grid_data;
kgrid.dt = 5e-8;
kgrid.t_array = [0:1500]*kgrid.dt;
input_args = {'Smooth', false, 'PMLInside', false, 'PlotPML', false};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% ******** Sending the reversed signals with detectors as time-verying sources ********
sensor_data1=fliplr(sensor_data);
nt1 = round(1.1*size(kgrid.t_array,2));
kgrid.t_array= [0:(nt1-1)]*kgrid.dt;
size(kgrid.t_array)
sensor.mask = 1*ones(Nx, Ny); % to store the temporal evolution of pressure at all grid points
sensor.mask(abs(Nx/2),abs(Ny/2))=1;
sensor.record = {'p'};
clear source
Fs = 1/kgrid.dt;
clear source
source.p_mask=grid_data;
source.p=[sensor_data1];
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);
P_data= sensor_data.p;
P_data_2D = reshape(P_data,Nx,Ny,[]);
figure,imagesc(((P_data_2D(:,:,1501)))), colormap hot, colorbar