I want to get some simulated data using the Kwave and reconstructed by the Time reversal method provided by the toolbox. However, some problem happened in the reconstructed image. There are some unexpected artifact appeared in the image. Could you please tell me why this happens? The following is the related code:
close all;
clear all;
clc;
addpath('C:\Program Files\MATLAB\R2013b\toolbox\k-Wave')
% =========================================================================
% SETUP THE GRID
% =========================================================================
% create the computational grid
Nx = 256; % number of grid points in the x (row) direction
Ny = 256; % number of grid points in the y (column) direction
dx = 0.04e-3; % grid point spacing in the x direction [m]
dy = dx; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);
center_freq = 5e6; % [Hz]
bandwidth = 80; % [%]
sensor.frequency_response = [center_freq, bandwidth];
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
% define the array of time points [s]
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
Nt = round(0.75*length(kgrid.t_array));
kgrid.t_array = (0:Nt-1)*dt;
% =========================================================================
% DEFINE A LINEAR ARRAY OF DIRECTIONAL ELEMENTS
% =========================================================================
array_left = zeros(Nx, Ny);
ut_left_x = 192;
i = 1;
for ut_left_y = 32:1:64
array_left(ut_left_x, ut_left_y) = 1;
ut_left_x = ut_left_x+1;
i = i+1;
end
sensor.mask = array_left;
figure, imagesc(sensor.mask); colormap(gray);
% =========================================================================
% DEFINE THE SOURCE
% =========================================================================
disc_magnitude = 10; % [au]
disc_radius = 1; % [grid points]
disc_sum = zeros(Nx, Ny);
disc_sum_x = zeros(Nx, Ny);
disc_sum_y = zeros(Nx, Ny);
for disc_y_pos = 64:30:192
disp(num2str(disc_y_pos));
for disc_x_pos = 64:30:192
disc_temp = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
disc_sum_x = disc_temp+ disc_sum_x;
end
disc_sum_y = disc_sum_x + disc_sum_y;
end
disc_sum = disc_sum_y;
disc_sum_index = find(disc_sum >= 1);
disc_sum(disc_sum_index) = 10;
figure, imagesc(disc_sum); colormap(gray);
source.p0 = disc_sum;
% smooth the initial pressure distribution and restore the magnitude
% source.p0 = smooth(kgrid, source.p0, true);
p0 = source.p0;
% =========================================================================
% RUN THE SIMULATION
% =========================================================================
PML_size = 30;
input_args = {'PMLInside', true, 'PMLSize', PML_size,'Smooth', false, 'PlotPML',false, 'PlotSim', true};
% run the simulation
tic
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
toc
%%%%%%%%%%%%%%%%%%%%% Time reversal reconstruction %%%%%%%%%%%%%%%%%%%%%%%%
% reset the initial pressure
source.p0 = 0;
% assign the time reversal data
sensor.time_reversal_boundary_data = sensor_data;
input_args = {'PMLInside', true, 'PMLSize', PML_size,'Smooth', false, 'PlotPML',false, 'PlotSim', true};
% run the time reversal reconstruction
p0_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% add first order compensation for only recording over a half plane
p0_recon = p0_recon*2;
% plot the initial pressure and sensor distribution
obj1 = p0;
obj1 = interp2(obj1,2);
% figure(1);imshow(obj1,[]);
figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, obj1);
colormap(hot); colorbar;
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
% plot the reconstructed initial pressure
img1 = double(p0_recon);
img1 = interp2(img1,2);
% figure(2);imshow(img1,[]);
figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, img1);
colormap(hot); colorbar;
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
img1_env = abs(hilbert(img1));
figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, img1_env);
colormap(hot); colorbar;
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
img1_log = 20*log10(img1_env);
img1_log = img1_log - max(img1_log(:));
figure; imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, img1_log, [-25,0]);
colormap(hot); colorbar;
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;