% 3D FFT Reconstruction For A Planar Sensor Example
%
% This example demonstrates the use of k-Wave for the reconstruction of a
% three-dimensional photoacoustic wave-field recorded over a planar sensor
% array. The sensor data is simulated using kspaceFirstOrder3D and
% reconstructed using kspacePlaneRecon. It builds on the Simulations In
% Three Dimensions and 2D FFT Reconstruction For A Line Sensor examples.
%
% author: Bradley Treeby
% date: 3rd 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;

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

% change scale to 2 to reproduce the higher resolution figures given in the
% examples
scale = 1;

% create the computational grid
PML_size = 8;       % size of the PML in pixels
Nx = 64*scale;      % number of pixels in the x direction
Ny = 64*scale;      % number of pixels in the y direction
Nz = 32*scale;      % number of pixels in the z direction
dx = 0.2e-3/scale;  % pixel width [m]
dy = 0.2e-3/scale;  % pixel length [m]
dz = 0.2e-3/scale;  % pixel height [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

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

% create initial pressure distribution using makeBall
ball_magnitude = 10;
ball_x_pos = 36*scale;    % pixels
ball_y_pos = 36*scale;    % pixels
ball_z_pos = 16*scale;    % pixels
ball_radius = 3*scale;    % pixels
source.p0 = ball_magnitude*makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

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

% define a binary planar sensor
sensor.mask = zeros(kgrid.Nz, kgrid.Nx, kgrid.Ny);
sensor.mask(1+PML_size, 1+PML_size:end-PML_size, 1+PML_size:end-PML_size) = 1;

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

% set the input arguements
input_args = {'PMLSize', PML_size, 'Smooth', false};

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

% reshape sensor data to x, y, t
[Np Nt] = size(sensor_data);
sensor_data_rs = reshape(sensor_data, sqrt(Np), sqrt(Np), Nt);

% reorder to t, x, y
sensor_data_rs = permute(sensor_data_rs, [3 1 2]);

% reconstruct the initial pressure
p_zxy = kspacePlaneRecon(sensor_data_rs, kgrid.dx, kgrid.dy, dt, medium.sound_speed, 'PosCond', true);

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

% trim the initial pressure distribution to remove the PML
source.p0 = source.p0(1+PML_size:end-PML_size, 1+PML_size:end-PML_size, 1+PML_size:end-PML_size);

% define slice locations
x_slice = ball_x_pos - PML_size;
y_slice = ball_y_pos - PML_size;
z_slice = ball_z_pos - PML_size;

% define a k-space grid using the dimensions of p_zxy
[Nz_recon Nx_recon Ny_recon] = size(p_zxy);
kgrid_recon = makeGrid(Nx_recon, kgrid.dx, Ny_recon, kgrid.dy, Nz_recon, dt*medium.sound_speed);

% define a k-space grid with the same z-spacing as the trimmed p0
[Nz_p0 Nx_p0 Ny_p0] = size(source.p0);
kgrid_interp = makeGrid(Nx_p0, kgrid.dx, Ny_p0, kgrid.dy, Nz_p0, kgrid.dz);

% resample the p_zxy to be the same size as p0
p_zxy_rs = interp3(kgrid_recon.x_off, kgrid_recon.z_off, kgrid_recon.y_off, p_zxy, kgrid_interp.x_off, kgrid_interp.z_off, kgrid_interp.y_off);

% plot the initial pressure
figure;
plot_scale = [-10 10];
subplot(2, 2, 1), imagesc(kgrid_interp.x(1,:,1)*1e3, kgrid_interp.z(:,1,1)*1e3, squeeze(source.p0(:, :, y_slice)), plot_scale);
title('x-z plane');
axis image;
subplot(2, 2, 2), imagesc(kgrid_interp.y(1,1,:)*1e3, kgrid_interp.z(:,1,1)*1e3, squeeze(source.p0(:, x_slice, :)), plot_scale);
title('y-z plane');
axis image;
xlabel('(All axes in mm)');
subplot(2, 2, 3), imagesc(kgrid_interp.x(1,:,1)*1e3, kgrid_interp.y(1,1,:)*1e3, squeeze(source.p0(z_slice, :, :)).', plot_scale);
title('x-y plane');
axis image;
colormap(getColorMap);

% plot the reconstructed initial pressure
figure;
subplot(2, 2, 1), imagesc(kgrid_interp.x(1,:,1)*1e3, kgrid_interp.z(:,1,1)*1e3, squeeze(p_zxy_rs(:, :, y_slice)), plot_scale);
title('x-z plane');
axis image;
subplot(2, 2, 2), imagesc(kgrid_interp.y(1,1,:)*1e3, kgrid_interp.z(:,1,1)*1e3, squeeze(p_zxy_rs(:, x_slice, :)), plot_scale);
title('y-z plane');
axis image;
xlabel('(All axes in mm)');
subplot(2, 2, 3), imagesc(kgrid_interp.x(1,:,1)*1e3, kgrid_interp.y(1,1,:)*1e3, squeeze(p_zxy_rs(z_slice, :, :)).', plot_scale);
title('x-y plane');
axis image;
colormap(getColorMap);

% view reconstruction slice by slice
flyThrough(p_zxy_rs);