% 3D Time Reversal Reconstruction For A Spherical Sensor Example
%
% This example demonstrates the use of k-Wave for the time-reversal
% reconstruction of a three-dimensional photoacoustic wave-field recorded
% over a spherical sensor. The sensor data is simulated and then
% time-reversed using kspaceFirstOrder3D. It builds on the 3D Time Reversal
% Reconstruction For A Planar Sensor Example. 
%
% author: Bradley Treeby
% date: 15th 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 = 10;       % size of the PML in pixels
Nx = 64;            % number of pixels in the x direction
Ny = 64;            % number of pixels in the y direction
Nz = 64;            % number of pixels in the z direction
dx = 0.2e-3;        % pixel width [m]
dy = 0.2e-3;        % pixel width [m]
dz = 0.2e-3;        % pixel height [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

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

% create initial pressure distribution using makeBall
ball_magnitude = 10;
ball_x_pos = 36;    % pixels
ball_y_pos = 32;    % pixels
ball_z_pos = 26;    % pixels
ball_radius = 3;    % pixels
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
p0 = smooth(p0, kgrid, true);

% define a Cartesian spherical sensor
sensor_radius = 4e-3;
num_sensor_points = 100;
center_pos = [0, 0, 0];
sensor_mask = makeCartSphere(sensor_radius, num_sensor_points, center_pos, true);

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

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

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

% assign the reconstruction grid to be the same
kgrid_recon = kgrid;

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

% run the time-reversal reconstruction
p0_recon = kspaceFirstOrder3D(sensor_data, kgrid_recon, c, rho, t_array, sensor_mask, input_args{:});

% create a binary sensor mask of an equivalent continuous sphere 
pixel_radius = round(sensor_radius/kgrid.dx);
binary_sensor_mask = makeSphere(kgrid_recon.Nx, kgrid_recon.Ny, kgrid_recon.Nz, pixel_radius);

% interpolate data to remove the gaps
[binary_sensor_data] = interpCartData(kgrid_recon, sensor_data, sensor_mask, binary_sensor_mask);

% run the time-reversal reconstruction
p0_recon_interp = kspaceFirstOrder3D(binary_sensor_data, kgrid, c, rho, t_array, binary_sensor_mask, input_args{:});

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

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

% plot the reconstructed initial pressure
figure;
plot_scale = [-0.2 0.2];
subplot(2, 2, 1), imagesc(kgrid.x(1,:,1)*1e3, kgrid.z(:,1,1)*1e3, squeeze(p0_recon(:, :, ball_y_pos)), plot_scale);
title('x-z plane');
axis image;
subplot(2, 2, 2), imagesc(kgrid.y(1,1,:)*1e3, kgrid.z(:,1,1)*1e3, squeeze(p0_recon(:, ball_x_pos, :)), plot_scale);
title('y-z plane');
axis image;
xlabel('(All axes in mm)');
subplot(2, 2, 3), imagesc(kgrid.x(1,:,1)*1e3, kgrid.y(1,1,:)*1e3, squeeze(p0_recon(ball_z_pos, :, :)).', plot_scale);
title('x-y plane');
axis image;
colormap(getColorMap);

% plot the reconstructed initial pressure
figure;
plot_scale = [-10 10];
subplot(2, 2, 1), imagesc(kgrid.x(1,:,1)*1e3, kgrid.z(:,1,1)*1e3, squeeze(p0_recon_interp(:, :, ball_y_pos)), plot_scale);
title('x-z plane');
axis image;
subplot(2, 2, 2), imagesc(kgrid.y(1,1,:)*1e3, kgrid.z(:,1,1)*1e3, squeeze(p0_recon_interp(:, ball_x_pos, :)), plot_scale);
title('y-z plane');
axis image;
xlabel('(All axes in mm)');
subplot(2, 2, 3), imagesc(kgrid.x(1,:,1)*1e3, kgrid.y(1,1,:)*1e3, squeeze(p0_recon_interp(ball_z_pos, :, :)).', plot_scale);
title('x-y plane');
axis image;
colormap(getColorMap);