% 2D Time Reversal Reconstruction For A Line Sensor Example
%
% This example demonstrates the use of k-Wave for the time-reversal
% reconstruction of a two-dimensional photoacoustic wave-field recorded
% over a linear sensor array. The sensor data is simulated and then
% time-reversed using kspaceFirstOrder2D. It builds on the 2D FFT
% Reconstruction For A Line Sensor Example. 
%
% author: Bradley Treeby
% date: 6th July 2009
% last update: 16th July 2009
%  
% This example is part of the k-Wave Toolbox (http://www.k-wave.org)

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
PML_size = 20;      % size of the PML in pixels
Nx = 256 - 2*PML_size; % number of pixels in the x (column) direction
Nz = 128 - 2*PML_size; % number of pixels in the z (row) direction
dx = 0.1e-3;        % pixel width [m]
dz = 0.1e-3;        % pixel height [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

% define the properties of the propagation medium
c = 1500;           % [m/s]
rho = 1000;         % [kg/m^3]

% create initial pressure distribution using makeDisc
disc_magnitude = 5;
disc_x_pos = 140;   % pixels
disc_z_pos = 60;    % pixels
disc_radius = 5;    % pixels
disc_2 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

disc_x_pos = 110;   % pixels
disc_z_pos = 30;    % pixels
disc_radius = 8;    % pixels
disc_1 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

p0 = disc_1 + disc_2;

% smooth the initial pressure distribution and restore the magnitude
p0 = smooth(p0, kgrid, true);

% define a binary line sensor
sensor_mask = zeros(kgrid.Nz, kgrid.Nx);
sensor_mask(1, :) = 1;

% create the time array
[t_array dt] = makeTime(kgrid, c);

% set the input arguements: force the PML to be outside the computational
% grid; switch off p0 smoothing within kspaceFirstOrder2D (as the grids for
% c and  rho in this example are homogeneous, the remaining smoothing
% parameters are not used so they can all be turned off)
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'Smooth', false};

% run the simulation
sensor_data = kspaceFirstOrder2D(p0, kgrid, c, rho, t_array, sensor_mask, input_args{:});

% add time reversal flag to the input options
input_args = [input_args, {'TimeRev', true}];

% run the time-reversal reconstruction
p0_recon = kspaceFirstOrder2D(sensor_data, kgrid, c, rho, t_array, sensor_mask, 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_zx = kspaceLineRecon(sensor_data.', kgrid.dx, dt, c, 'PosCond', true, 'Interp', '*cubic');

% define a second k-space grid using the dimensions of p_zx
[Nz_recon Nx_recon] = size(p_zx);
kgrid_recon = makeGrid(Nx_recon, kgrid.dx, Nz_recon, dt*c);

% resample the p_zx to be the same size as p0
p_zx_rs = interp2(kgrid_recon.x_off, kgrid_recon.z_off, p_zx, kgrid.x_off, kgrid.z_off);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0 + sensor_mask*disc_magnitude, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;
colorbar;

% plot the reconstructed initial pressure 
figure;
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0_recon, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-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.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0_recon, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;
colorbar;

% plot a profile for comparison
figure;
plot(kgrid.x(1,:)*1e3, p0(disc_z_pos, :), 'k-', kgrid.x(1,:)*1e3, p_zx_rs(disc_z_pos, :), 'r--',  kgrid.x(1,:)*1e3, p0_recon(disc_z_pos, :), 'b:');
xlabel('x-position [mm]');
ylabel('Pressure');
legend('Initial Pressure', 'FFT Reconstruction', 'Time Reversal');
axis tight;
set(gca, 'YLim', [0 5.1]);