I modified "2D Time Reversal Reconstruction For A Line Sensor Example" and made medium heterogeneous with the following addition: However, simulation results in an expected situation. Did I do something wrong? Is it possible to apply image reconstruction for defects?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(I added an air bubble into a steel block)
disc_x_pos=round(Nx*0.7);
disc_y_pos=round(Ny/4);
disc_radius=5;
air_buble = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
air_buble_reverse = (air_buble-1)*-1;
medium.sound_speed = medium.sound_speed .* air_buble_reverse;
medium.sound_speed = medium.sound_speed + air_buble * 340 ;
medium.density = 7800*ones(Nx, Ny); % [kg/m^3]
medium.density = medium.density .* air_buble_reverse;
medium.density = medium.density + air_buble * 12.25 ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
therefore, the complete solution is:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
% =========================================================================
% 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);
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [au]
disc_x_pos = 60; % [grid points]
disc_y_pos = 140; % [grid points]
disc_radius = 5; % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
disc_x_pos = 30; % [grid points]
disc_y_pos = 110; % [grid points]
disc_radius = 8; % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
medium.sound_speed = 1500*ones(Nx, Ny); % [m/s]
disc_x_pos=round(Nx*0.7);
disc_y_pos=round(Ny/4);
disc_radius=5;
air_buble = makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);
air_buble_reverse = (air_buble-1)*-1;
medium.sound_speed = medium.sound_speed .* air_buble_reverse;
medium.sound_speed = medium.sound_speed + air_buble * 340 ;
medium.density = 7800*ones(Nx, Ny); % [kg/m^3]
medium.density = medium.density .* air_buble_reverse;
medium.density = medium.density + air_buble * 12.25 ;
% smooth the initial pressure distribution and restore the magnitude
%p0 = smooth(kgrid, disc_1 + disc_2, true);
p0 = smooth(kgrid, disc_1 + disc_2, true);
% assign to the source structure
source.p0 = p0;
% define a binary line sensor
sensor.mask = zeros(Nx, Ny);
sensor.mask(1, :) = 1;
% create the time array
[kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);
% 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};
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% reset the initial pressure
source.p0 = 0;
% assign the time reversal data
sensor.time_reversal_boundary_data = sensor_data;
% 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;
% repeat the FFT reconstruction for comparison
p_xy = kspaceLineRecon(sensor_data.', dy, dt, medium.sound_speed, 'PosCond', true, 'Interp', 'linear');
% define a second k-space grid using the dimensions of p_xy
[Nx_recon, Ny_recon] = size(p_xy);
kgrid_recon = makeGrid(Nx_recon, dt*medium.sound_speed, Ny_recon, dy);
% resample p_xy to be the same size as source.p0
p_xy_rs = interp2(kgrid_recon.y, kgrid_recon.x - min(kgrid_recon.x(:)), p_xy, kgrid.y, kgrid.x - min(kgrid.x(:)));
% =========================================================================
% VISUALISATION
% =========================================================================
% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0 + sensor.mask*disc_magnitude, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;
% plot the reconstructed initial pressure
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0_recon, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;
% apply a positivity condition
p0_recon(p0_recon < 0) = 0;
% plot the reconstructed initial pressure with positivity condition
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0_recon, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;
% plot a profile for comparison
figure;
plot(kgrid.y_vec*1e3, p0(disc_x_pos, :), 'k-', ...
kgrid.y_vec*1e3, p_xy_rs(disc_x_pos, :), 'r--', ...
kgrid.y_vec*1e3, p0_recon(disc_x_pos, :), 'b:');
xlabel('y-position [mm]');
ylabel('Pressure');
legend('Initial Pressure', 'FFT Reconstruction', 'Time Reversal');
axis tight;
set(gca, 'YLim', [0 5.1]);