Dear Bradley:
Thank you for your help. I am trying to reconstruct a disc medium density model with two steps: 1. forward modeling (a. with disc density model, b. without disc density model, c. data is a-b), and 2. time reversal to reconstruct the medium density model. I tried to use a pulse source with the Mexican hat wavelet. I modified the 2D Time Reversal Reconstruction For A Line Sensor Example as below. I keep all the timing parameters same in forward and reversal modeling. However, the reconstructed image is not same as the disc density model (it seems reconstructed time is a little longer than it should be). I was wondering if you could please help with my questions?
My questions are as follows:
1. The time parameters in reconstruction are kept same as those in the forward modeling, could you please help with the possible reason why the reconstructed image is not the same as the medium disc2? How to restore the exact the original medium disc2 shape (e.g. do I need to control the timing)?
2. the reconstructed image 'density_recon' values have both positive and negative values,
why there are negative image values (the medium disc2 are postive values as medium density? Why the original disc values can't be restored?
"density_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % inverted medium parameters"
Thank you very much for your help in advance and best regards,
John
p.s. the Matlab script is as below:
--------------------------------------------------------------------------------
clear all;
% =========================================================================
% SET UP THE SIMULATION
% =========================================================================
% create the computational grid
PML_size = 20; % size of the PML in grid points
Nx = 128 - 2*PML_size; % number of grid points in the x (row) direction
Ny = 256 - 2*PML_size; % number of grid points in the y (column) direction
dx = 0.1e-3; % grid point spacing in the x direction [m]
dy = 0.1e-3; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, 1500);
% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [au]
disc_x_pos = 10; % [grid points]
disc_y_pos = floor(Ny/2); % [grid points]
disc_radius = 5; % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
% smooth the initial pressure distribution and restore the magnitude
% scale the initial pressure distribution
p0_magnitude = 1;
p0 = smooth(kgrid, disc_1, true);
% assign to the source structure
source_freq = 10;
source.p_mask = disc_1;
source.p = p0_magnitude .* (1 - 2 * pi^2 * source_freq^2 .* (kgrid.t_array).^2)...
.* exp(-pi^2 * source_freq^2 .* (kgrid.t_array).^2); % Ricker wavelet
% define the properties of the propagation medium
disc_x_pos = floor(Nx/2); % [grid points]
disc_y_pos = floor(Ny/2); % [grid points]
disc_radius = 5; % [grid points]
disc_2 = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
p3 = smooth(kgrid, disc_2, true);
density_input = 1000*ones(Nx, Ny); % [kg/m^3]
medium.sound_speed = 1500 ; % [m/s]
medium.density = density_input + 1000 * disc_2; % [m/s]
% set the input arguements: force the PML to be outside the computational
% grid; switch off p0 smoothing within kspaceFirstOrder2D
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false, 'PlotPML', false};
% =========================================================================
% DEFINE SENSOR INDICES FOR LATER USE
% =========================================================================
% define a four-sided box array and find the indices of the sensors
sensor.mask = zeros(Nx, Ny);
sensor.mask(1,:) = 1;
sensor.mask(:,1) = 1;
sensor.mask(end,:) = 1;
sensor.mask(:,end) = 1;
sensor_indices = find(sensor.mask==1);
% find the indices along two sides, with respect to sensor_indices
sensor.mask = zeros(Nx, Ny);
sensor.mask(end,:) = 1;
sensor.mask(:,end) = 1;
sensor_box_indices1 = find(sensor.mask(sensor_indices)==1);
% define an L-shaped array used for the measurements, find the sensor
% indices, and their indices with respect to sensor_indices
sensor.mask = zeros(Nx, Ny);
sensor.mask(1,:) = 1;
sensor.mask(:,1) = 1;
original_sensor_indices = find(sensor.mask == 1);
sensor_box_indices2 = find(sensor.mask(sensor_indices)==1);
% =========================================================================
% RUN THE SIMULATION AND RECORD THE DATA OVER THE L-SHAPED ARRAY
% =========================================================================
sensor_data_object = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
medium.density = density_input; % [m/s]
sensor_data_no_object = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
sensor_data = sensor_data_object - sensor_data_no_object;
% =========================================================================
% RECONSTRUCT AN IMAGE USING TIME REVERSAL AND LINE ARRAY DATA
% =========================================================================
% define a line array of sensors
sensor.mask = zeros(Nx, Ny);
sensor.mask(1, :) = 1;
% assign the time reversal data
sensor_data = sensor_data; % - sensor_data_no_object;
sensor.time_reversal_boundary_data = 1 * sensor_data((sensor.mask(original_sensor_indices) == 1),:);
% set the initial pressure to be zero
source.p = source.p * 0;
% run the time reversal reconstruction
density_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % inverted medium parameters
figure;
imagesc(disc_2)
title('medium density');
colorbar;
figure;
imagesc(sensor_data)
title('sensor_data');
colorbar;
figure;
imagesc(density_recon)
title('density_recon');
colorbar;