% Image Reconstruction With Bandlimited Sensors Example
%
% This example demonstrates how the bandwidth of sensor elements can give
% rise to artefacts in time reversal photoacoustic image reconstruction. It
% builds on the previous 2D time reversal examples. 
%
% author: Bradley Treeby
% date: 27th August 2010
% last update: 9th September 2010
%  
% 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
Nx = 128;           % number of pixels in the x (column) direction
Nz = 128;           % 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
medium.sound_speed = 1500;  % [m/s]

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

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [au]
disc_x_pos = 50;    % [pixels]
disc_z_pos = 50;    % [pixels]
disc_radius = 8;    % [pixels]
disc_1 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

disc_magnitude = 4; % [au]
disc_x_pos = 85;    % [pixels]
disc_z_pos = 65;    % [pixels]
disc_radius = 7;    % [pixels]
disc_2 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

disc_magnitude = 3; % [au]
disc_x_pos = 60;    % [pixels]
disc_z_pos = 80;    % [pixels]
disc_radius = 5;    % [pixels]
disc_3 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

source.p0 = disc_1 + disc_2 + disc_3;

% define a centered circular sensor
sensor_radius = 44;   % [pixels]
sensor.mask = makeCircle(Nx, Nz, Nx/2, Nz/2, sensor_radius);

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% filter the sensor data using a high pass filter
Fs = 1/dt;              % [Hz]
cutoff_freq = 1e6;      % [Hz]
sensor_data_low_pass = zeros(size(sensor_data));
for index = 1:sum(sensor.mask(:))
    sensor_data_low_pass(index, :) = applyFilter(sensor_data(index, :), Fs, cutoff_freq, 'HighPass', 'ZeroPhase', true);
end

% filter the sensor data using a low pass filter
Fs = 1/dt;              % [Hz]
cutoff_freq = 1e6;      % [Hz]
sensor_data_high_pass = zeros(size(sensor_data));
for index = 1:sum(sensor.mask(:))
    sensor_data_high_pass(index, :) = applyFilter(sensor_data(index, :), Fs, cutoff_freq, 'LowPass', 'ZeroPhase', true);
end

% assign the sensor frequency response
center_freq = 3e6;      % [Hz]
bandwidth = 100;        % [%]
sensor.frequency_response = [center_freq, bandwidth];

% re-run the simulation to obtain the filtered sensor data
sensor_data_gaussian = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% =========================================================================
% TIME REVERSAL IMAGE RECONSTRUCTION
% =========================================================================

% reset the initial pressure
source.p0 = 0;

% assign the unfiltered time reversal data
sensor.time_reversal_boundary_data = sensor_data;

% run the time reversal reconstruction
p0_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% assign the high pass filtered time reversal data
sensor.time_reversal_boundary_data = sensor_data_low_pass;

% re-run the time reversal reconstruction
p0_recon_low_pass = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% assign the low pass filtered time reversal data
sensor.time_reversal_boundary_data = sensor_data_high_pass;

% re-run the time reversal reconstruction
p0_recon_high_pass = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% assign the gaussian filtered time reversal data
sensor.time_reversal_boundary_data = sensor_data_gaussian;

% re-run the time reversal reconstruction
p0_recon_gaussian = kspaceFirstOrder2D(kgrid, medium, source, sensor);

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

% plot the reconstructed initial pressure with no filtering
figure;
mx = max(abs(p0_recon(:)));
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0_recon + mx*sensor.mask, [-mx, mx]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot the reconstructed initial pressure with high pass filtering
figure;
mx = max(abs(p0_recon_low_pass(:)));
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0_recon_low_pass + mx*sensor.mask, [-mx, mx]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot the reconstructed initial pressure with low pass filtering
figure;
mx = max(abs(p0_recon_high_pass(:)));
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0_recon_high_pass + mx*sensor.mask, [-mx, mx]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot the reconstructed initial pressure with Gaussian filtering
figure;
mx = max(abs(p0_recon_gaussian(:)));
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0_recon_gaussian + mx*sensor.mask, [-mx, mx]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

% plot a profile through the two reconstructions
figure;
x = kgrid.x(1,:)*1e3;
plot(x, p0_recon(50, :), 'k-', x, p0_recon_low_pass(50, :), 'b-', x, p0_recon_high_pass(50, :), 'g-', x, p0_recon_gaussian(50, :), 'r-');
set(gca, 'XLim', [min(x), max(x)]);
ylabel('Pressure [au]');
xlabel('x-position [mm]');
legend('Not Filtered', 'High Pass Filtered', 'Low Pass Filtered', 'Gaussian Filtered');
