% 2D Time Reversal Reconstruction For A Circular 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 circular sensor array. The sensor data is simulated and then
% time-reversed using kspaceFirstOrder2D. It builds on the 2D Time Reversal
% Reconstruction For A Line Sensor Example. 
%
% author: Bradley Treeby
% date: 7th July 2009
% last update: 19th January 2010
%  
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>. 

clear all;

% load the initial pressure distribution from an image and scale
p0_magnitude = 2;
p0 = p0_magnitude*loadImage('EXAMPLE_source_two.bmp');

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

% resize the input image to the desired number of pixels
p0 = resize(p0, Nx, Nz);

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

% assign to the source structure
source.p0 = p0;

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.density = 1000;   	% [kg/m^3]

% define a centered Cartesian circular sensor
sensor_radius = 4.5e-3;     % [m]
sensor_angle = 3*pi/2;      % [rad]
sensor_pos = [0, 0];
num_sensor_points = 70;
cart_sensor_mask = makeCartCircle(sensor_radius, num_sensor_points, sensor_pos, sensor_angle);

% assign to sensor structure
sensor.mask = cart_sensor_mask;

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

% set the input options
input_args = {'Smooth', [false, true], 'PMLInside', false};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% add noise to the recorded sensor data
sensor_data = sensor_data + max(sensor_data(:))*0.05*(rand(size(sensor_data)) - 0.5);

% create a second computation grid for the reconstruction to avoid the
% inverse crime
Nx = 300;           % number of pixels in the x (column) direction
Nz = 300;           % number of pixels in the z (row) direction
dx = x/Nx;          % pixel width [m]
dz = z/Nz;          % pixel height [m]
kgrid_recon = makeGrid(Nx, dx, Nz, dz);

% attach the original time array
kgrid_recon.t_array = kgrid.t_array;

% 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_recon, medium, source, sensor, input_args{:});

% create a binary sensor mask of an equivalent continuous circle 
pixel_radius = round(sensor_radius/kgrid_recon.dx);
binary_sensor_mask = makeCircle(kgrid_recon.Nx, kgrid_recon.Nz, kgrid_recon.Nx/2, kgrid_recon.Nz/2, pixel_radius, sensor_angle);

% assign to sensor structure
sensor.mask = binary_sensor_mask;

% interpolate data to remove the gaps and assign to sensor structure
sensor.time_reversal_boundary_data = interpCartData(kgrid_recon, sensor_data, cart_sensor_mask, binary_sensor_mask);

% run the time-reversal reconstruction
p0_recon_interp = kspaceFirstOrder2D(kgrid_recon, medium, source, sensor, input_args{:});

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

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

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;

% plot the reconstructed initial pressure 
figure;
imagesc(kgrid_recon.x(1,:)*1e3, kgrid_recon.z(:,1)*1e3, p0_recon, [-1 1]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot the reconstructed initial pressure using the interpolated data
figure;
imagesc(kgrid_recon.x(1,:)*1e3, kgrid_recon.z(:,1)*1e3, p0_recon_interp, [-1 1]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot a profile for comparison
slice_pos = 4.5e-3;  % [m] location of the slice from top of grid [m]
figure;
plot(kgrid.x(1,:)*1e3, p0(round(slice_pos/kgrid.dx), :), 'k--', kgrid_recon.x(1,:)*1e3, p0_recon(round(slice_pos/kgrid_recon.dx), :), 'r-',  kgrid_recon.x(1,:)*1e3, p0_recon_interp(round(slice_pos/kgrid_recon.dx), :), 'b-');
xlabel('x-position [mm]');
ylabel('Pressure');
legend('Initial Pressure', 'Point Reconstruction', 'Interpolated Reconstruction');
axis tight;
set(gca, 'YLim', [0 2.1]);