% 2D FFT Reconstruction For A Line Sensor Example
%
% This example demonstrates the use of k-Wave for the reconstruction of a
% two-dimensional photoacoustic wave-field recorded  over a linear array of
% sensor elements  The sensor data is simulated using kspaceFirstOrder2D
% and reconstructed using kspaceLineRecon. It builds on the Homogeneous
% Propagation Medium and Heterogeneous Propagation Medium examples. 
%
% author: Bradley Treeby
% date: 2nd July 2009
% last update: 7th February 2011
%  
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010, 2011 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;

% =========================================================================
% 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 direction
Nz = 128 - 2*PML_size;  % number of pixels in the z 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
medium.sound_speed = 1500;           % [m/s]

% 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);

source.p0 = disc_1 + disc_2;

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

% define a binary line sensor
sensor.mask = zeros(kgrid.Nz, kgrid.Nx);
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, 'PlotPML', false, 'Smooth', false};

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

% reconstruct the initial pressure
p_zx = kspaceLineRecon(sensor_data.', kgrid.dx, dt, medium.sound_speed, 'PlotRecon', true, 'PosCond', true);

% 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*medium.sound_speed);

% resample p_zx to be the same size as source.p0
p_zx_rs = interp2(kgrid_recon.x, kgrid_recon.z - min(kgrid_recon.z(:)), p_zx, kgrid.x, kgrid.z - min(kgrid.z(:)));

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

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

% 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.x(1,:)*1e3, kgrid.z(:,1)*1e3, p_zx_rs, [-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, source.p0(disc_z_pos, :), 'k-', kgrid.x(1,:)*1e3, p_zx_rs(disc_z_pos, :), 'r--');
xlabel('x-position [mm]');
ylabel('Pressure');
legend('Initial Pressure', 'Reconstructed Pressure');
axis tight;
set(gca, 'YLim', [0 5.1]);