% Quantitative Microwave Holography 
% Written by Daniel Tajik, McMaster University
% Written to complement:
%   Introduction to Microwave Imaging, by Natalia K. Nikolva
%   Section 5.2: 3D holographic inversion with wideband planar data acq.
%
%   Example included is based on section 5.7
% 
% Current Revision 1 (July 5th, 2017)
% Changes:

% Clear all previous information every execution, and calculate runtime.
clc;
clear all;
close all;
tic

%% Initial Information: Loading data files, declaring variables
%--------------------------------------------------------------------------
%---------------------------USER INPUT BELOW-------------------------------
%-------------------------LOAD YOUR DATA HERE------------------------------
%--------------------------------------------------------------------------
%-------------------------SIMULATION EXAMPLE-------------------------------
%--------------------------------------------------------------------------

% EXAMPLE
% This code can be run with provided S-parameter data somulated in the
% F-shape example; see Section 4.3.5.

% Begin loading the files containing the RO (reference object), CO
% (calibration object) and OUT (object-under-test) datasets. All datasets
% are acquired using S-parameters
% Here, S11 = Reflection coefficient at port 1
%       S21 = Transmission coefficient from port 1 to port 2
%       S12 = Transmission coefficient from port 2 to port 1
%       S22 = Reflection coefficient at port 2



% LOAD AND ORGANIZE DATASETS
%--------------------------------------------------------------------------
% RO = Reference object (only the background medium; no embedded scatterers)
load S11_RO.mat
load S21_RO.mat
load S12_RO.mat
load S22_RO.mat
 
% Calibration object (background medium with a single scatterer; the point-
% spread function)
load S11_CO.mat
load S21_CO.mat
load S12_CO.mat
load S22_CO.mat
 
% Object under test
load S11_OUT.mat
load S21_OUT.mat
load S12_OUT.mat 
load S22_OUT.mat

% To simplify the code further on, the S-parameters are placed within cell
% structures. The order between the RO, CO and OUT MUST be aligned with the
% appropriate reflection/transmission coefficient

s_RO = {S11_RO;S21_RO;S12_RO;S22_RO};
s_CO = {S11_CO;S21_CO;S12_CO;S22_CO};
s_OUT = {S11_OUT;S21_OUT;S12_OUT;S22_OUT};


% DECLARE VARIABLES FOR QUANTITATIVE RECONSTRUCTION
%--------------------------------------------------------------------------
% (UNIT OF MEASUREMENT = METERS)

% Relative permittivity of background medium
eps_RO = 1;

% Relative permittivity of scattering probe in CO
eps_CO = 1.1;

% Volume of the scattering probe 
vol_sp = 10e-3^3; 
% CO_area = 5e-3^2;

% Scan step
xstep = 10e-3;

% Range of CO measurements along z. In the case of 2D imaging, just put the
% distance between the antennas
z = 30e-3:10e-3:50e-3;

% Bandwidth of frequencies
ftable = 3e9:1e9:16e9;



% ADDITIONAL SETTINGS 
%--------------------------------------------------------------------------
% Replace 0 with 1 to enable the boolean variables

% The default approximation used in this code is Born's approximation (see
% section 1.14). However, this method can also operate using Rytov's
% approximation (see section 1.16). Setting the parameter to 1 will execute
% the code using Rytov's approximation.
enable_rytov = 0;

% NOT REQUIRED FOR SIM EXAMPLE
% Apodization filtering: This filter will alleviate errors related to the
% DFT as shown in the work of [Tajik, 2017, submitted]. While it has no
% impact on simulated data, it is typically required when processing
% experimental data.
apod = 0;
% Shows the apodization filter.
show_apod = 0;

% NOT REQUIRED FOR SIM EXAMPLE
% Post-inversion low-pass filtering: This filter will alleviate errors
% related to oversampling (see section 5.3). In the simulation example 
% provided, this is not the case, and is not necessary. In experimental 
% data, this filter is a necessity.
lpf = 0;
% Shows the apodization filter.
show_lpf = 0;

%Smoothing filter: This increases the sampling rate of the final image, and
%smooths it with a low pass filter. Use if you would like to see a smoother
%result.
smoothing = 0;

% Save the images if desired
save_comb = 0;
save_ind = 0;
% The images can either be generated per each S-parameter, or all together.
% The simulation requires all 4 S-parameters.
close_comb = 0;
close_ind = 1;

% PARAMETER DELCARATIONS 
%--------------------------------------------------------------------------
% Sets appropriate parameters used later on in the code. No adjustment 
% necessary unless specifically required.

% Removes the background RO from the CO and OUT data. If Rytov is selected,
% it will perform the Rytov version of de-embedding (see [Tajik, IMS 2017])
subtract_RO = [1;1;1;1];

% declaring required parameters used further on in the code.
ns = size(s_CO,1); 
[nx ny nf nz] = size(s_CO{1});
c0 = 2.998e8;
xrng = xstep*(nx-1);
x=linspace(-xrng/2,xrng/2,nx);
y=linspace(-xrng/2,xrng/2,ny);
% if it is a 2D image, zstep is the distance between antennas.
if length(z)>1
    zstep = z(2)-z(1);
else
    zstep = z(1);
end
units = 'm';
 
%--------------------------------------------------------------------------
%---------------------------END OF USER INPUT------------------------------
%------------TREAD CAUTIOUSLY WHEN EDITTING BEYOND THIS POINT--------------
%--------------------------------------------------------------------------

%% Stage One: Quantitative Microwave Holography
%--------------------------------------------------------------------------
%The basics of this algorithm can be seen in Section 5.2 of Introduction to 
%Microwave Imaging. It has the following stages:
%1.1 Acquisition of scattering parameter data
%1.2 Apodization Filtering
%1.3 Normalization and the Fast Fourier Transform 
%1.4 Inversion step at each spectral position k
%1.5 Post inversion low pass filtering
%1.6 Multiplicative factoring and final image adjustments
%1.7 Image Generation
%--------------------------------------------------------------------------

%% Part 1.1: Acquiring the scattered response (Born or Rytov Approximation)
% 
% The scattered response is derived from the CO and OUT responses

for i=1:length(s_RO) 
    if size(s_RO{i},1) == 1
        for f = 1:nf
            % in our simulation data, the RO is typically one measurement
            % and needs to be built into a map, which is not the case for
            % experimental measurements
            s_RO_hold{i}(:,:,f) = s_RO{i}(1,1,f)*ones(nx,ny);
        end
        s_RO{i} = s_RO_hold{i};
    end
    
    
    %BORN STRATEGY
    %----------------------------------------------------------------------
    % Scattered response is S_SC = S_(OUT,CO) - S_RO
    if subtract_RO(i) == 1 && ~enable_rytov
        %loop through each frequency
        for f = 1:nf
            s_OUT{i}(:,:,f) = s_OUT{i}(:,:,f) - s_RO{i}(:,:,f,1);

            %loop through each z-plane
            for iz = 1:nz
                s_CO{i}(:,:,f,iz) = s_CO{i}(:,:,f,iz) - s_RO{i}(:,:,f);
            end
        end
    end
end

    %RYTOV STRATEGY
    %----------------------------------------------------------------------
    % Scattered response is S_SC = S_RO*ln(S_(OUT,CO)/S_RO)
if enable_rytov==1
    if isempty(s_RO)
        disp('No RO was loaded! Output will be incorrect.')
        enable_rytov = 0;
    else
        disp('Solving using the Rytov Approximation')
        %loop through each s parameter
        for n = 1:ns
            s_OUT{n} = s_RO{n}.*(log(s_OUT{n}./s_RO{n}));
 
            %loop through each z plane
            for iz = 1:nz 
                s_CO{n}(:,:,:,iz) =  s_RO{n}(:,:,:,1).*(log(squeeze(s_CO{n}(:,:,:,iz))./s_RO{n}(:,:,:,1)));
            end          
        end
    end
else
    disp('Solving using the Linear Born Approximation')
end

%% Part 1.2: Apodization Filtering
%--------------------------------------------------------------------------
% Important if there are large amounts of noise/discontinuities at the
% edges of the ROI. Typically left disabled in simulation results. Several
% filters are included, but using anything beyond the cosine filter
% requires specific adjustments.
%--------------------------------------------------------------------------
fftw('planner','exhaustive');

if  apod==1
    
    %cosine filter. This is the default apodization filter
    %----------------------------------------------------------------------
    w1 = cos(linspace(-pi/2, pi/2, nx));
    w2 = cos(linspace(-pi/2, pi/2, ny));
    cos_filt = w1'*w2;
 
    
   %other filters
   %-----------------------------------------------------------------------
    %There are two extra filters that you may wish to try, a gaussian
    %filter and the Happ-Genzel filter. They may need to be manually
    %adjusted in order to function corectly.
%    apod_adj = 1;
%    adj = apod_adj;
%    stddev = nx/2*apod_adj;
%     apod_filter = zeros(nx);
%     for ix = 1:nx
%         for iy = 1:ny
%             
%             %gaussian filter
%             g_filt(ix,iy) = exp((-(ix-nx/2).^2)/(2*(stddev/2)^2))*exp((-(iy-ny/2).^2)/(2*(stddev/2)^2));
%             
%             %Happ-Genzel filter
%             happ_genzel(ix,iy) = (0.54+0.46*cos(pi*(ix+nx/2)/(nx/2)))*(0.54+0.46*cos(pi*(iy+ny/2)/(ny/2)));
%         end    
%     end
%     apod_filter = g_filt;
%     apod_filter = happ_genzel;
%   %-----------------------------------------------------------------------


apod_filter = cos_filt;
else
    %if we dont want a gauss filter, just multiply everything by 1!
    apod_filter = ones(nx);
end

 %show the apodization filter
if show_apod
    figure;
    imagesc(apod_filter);
    title('Apodization Filter');
end

%% Part 1.3: Normalization and the Fast Fourier Transform 
%--------------------------------------------------------------------------
% Compute the fourier transforms for the CO and OUT. We also normalize the
% signals to avoid biasing certain frequency values. The apodization is
% applied before the FFT
%--------------------------------------------------------------------------

%FFT is computed independently on each (x,y) image. 
% loop through each S-parameter
for n = 1:ns
    SF_CO{n} = zeros(nx,ny,nf,nz);
    SF_OUT{n} = zeros(nx,ny,nf);

    %loop through each frequency
    for fi = 1:nf
        %calculate the max value in the CO measurement at each frequency to
        %use in normalization
        maxi = max(max(max(abs(s_CO{n}(:,:,fi,:)))));
        
        %OUT FFT and normalization
        SF_OUT{n}(:,:,fi)=fft2(s_OUT{n}(:,:,fi).*apod_filter/maxi); 
           
        %CO FFT and normalization looped through each z plane
        for iz = 1:nz            
            SF_CO{n}(:,:,fi,iz)=fft2((s_CO{n}(:,:,fi,iz).*apod_filter/maxi));
        end
        
    end
   
end

%% Part 1.4: Solving the Inverse Problem: Least Squares Solution
%--------------------------------------------------------------------------
% Solve the QMH least squares problem. We go through each voxel and solve
% the least squares problem independently across the combinations of
% S-parameters, as well as combining the S-parameters together to find the
% minimum 2 norm solution. See equation 5.21
%
% Independet measurements are indentified by the suffix _ind
% Combined min 2 norm results are identified by the suffix _tot
% 
% Be careful when looking at indepedent data... it can generate strange
% values in simulated datasets! Especially when you have linear dependent
% CO planes
%--------------------------------------------------------------------------

% at every position (kx,ky) the inversion is performed
for kx = 1:nx
    for ky = 1:ny
        
        A_ind = [];
        b_ind = [];
        A_tot = [];
        b_tot = [];

        %This loop process each S-parameter independently
        for n = 1:ns

            A_ind = reshape(SF_CO{n}(kx,ky,:,:),[nf,nz]);
            b_ind = reshape(SF_OUT{n}(kx,ky,:),[nf, 1]);
            
            %The S-parameters are stacked to process the result processed
            %using all captured data. 
            
            A_tot = [A_tot;A_ind];
            b_tot = [b_tot;b_ind];
 
            %Perform inversion on ONLY 1 S-parameter
            Fout_ind(kx,ky,:,n) = pinv(A_ind)*(b_ind);

        end

        %Perform inversion on the combined S-parameters
        Fout_tot(kx,ky,:) = pinv(A_tot)*b_tot;   
    end
end

%% Part 1.5: Post-inversion low pass filtering
%--------------------------------------------------------------------------
% Filtering has been show to be a vital part of the process of quantitative
% and qualititative holography, especially when there is oversampling
% involved. Now we will generate a Butterworth filter to remove high
% frequency noise reponses. See [Tajik, 2017, submitted] for more
% information.
% The filter sets the limit of the fourier space equal to the limit as
% stated in section 5.3 of Introduction of Microwave Imaging
%--------------------------------------------------------------------------

if lpf == 1
   
    %calculate minimum wavelength
    lambda_min = c0/sqrt(abs((eps_RO)))/ftable(nf);
    
    %Calculate the maximum angle the antenna can see. This assumes that the
    %directivity of the antenna does not affect the maximum angle, which is
    %typically not the case in experimental reconstructions
    alpha = atan(sqrt(2*(xrng)^2)/z(nz));
    
    %calculate maximum kmax, and determine is location in the physical
    %FFT images
    kmax = 4*pi*sin(alpha)/lambda_min;
    kx = pi/xstep;
    dkx = 2*kx/nx;
    kmax_steps = kmax/dkx;
    dev =kmax_steps;
    
    %set the orde of the Butterworth filter
    order=2;
    bfilt = ones(nx,ny);
    %BUTTERWORTH FILTER
    for i = 1:nx
        for j =1:ny
            chg = ((i-((nx-1)/2+1))^2 + (j-((ny-1)/2+1))^2)^.5;
            % Create Butterworth filter.
            bfilt(i,j)= 1/(1 + (chg/dev)^(2*order));
        end
    end
    %Shift the filter to align with the FFT versions of the CO and OUT
    lpf = fftshift(bfilt);
else
    %if we dont want to filter, just multiply everything by 1!
    lpf = ones(nx);
end


 %show the filter
if show_lpf
    figure;
    imagesc(ifftshift(lpf));
    title('Post-Inv. L.P. Filter');
end

%% Part 1.6: Multiplicative factoring and final image adjustments
%--------------------------------------------------------------------------
%This is the multiplication factor that produces the quantitative
%information. Again, see section 5.3
%--------------------------------------------------------------------------

%this factor mulitplies with the reconstruction results. See equation 5.22
chg_eps = (vol_sp/(xstep^2*zstep))*(eps_CO-eps_RO);

for iz = 1:nz

    %smoothing for the combined S-parameter data.
    if smoothing == 1

        %ratio at which it is resampled
        p = 4;
        q = 1;  
   
        %has to be done one dimension at a time
        first_dim = resample(fftshift(ifft2(lpf.*Fout_tot(:,:,iz))),p,q);
        both_dim = resample(first_dim',p,q)';
        h = 1/4*[1 1;1 1];

        fout_tot(:,:,iz) = imfilter(both_dim,h);
    else
        %Typically, no smoothing is applied
        fout_tot(:,:,iz) = (fftshift(ifft2(lpf.*Fout_tot(:,:,iz))));
    end
    
    % Apply the multiplicative processing to obtain quantitative
    % information. See equation 5.22
    epsout_tot(:,:,iz) = chg_eps*fout_tot(:,:,iz) + eps_RO;
    
    
    % The above is repeated for each independent S-parameter
    for n = 1:ns
        

        if smoothing == 1
            
            first_dim = resample(fftshift(ifft2(lpf.*(Fout_ind(:,:,iz,n)))),p,q);
            both_dim = resample(first_dim',p,q)';
            
            h = 1/4*[1 1;1 1];
            
            fout_ind(:,:,iz,n) = imfilter(both_dim,h);
        else
            %typically no smoothing is required! Apply the 
            fout_ind(:,:,iz,n) = (fftshift(ifft2(lpf.*(Fout_ind(:,:,iz,n)))));
        end       
        
        %equation 5.22
        epsout_ind(:,:,iz,n) = chg_eps*fout_ind(:,:,iz,n) + eps_RO;
        end
    
end

%% Part 1.7: Image Generation
%Final images are generated here.


for iz = 1:nz
    for n = 1:ns

        
        maintitle = sprintf('S-param = %d, Z-plane = %.02f %s',n,z(iz),units);
        filename = sprintf('S%d_z%.03f.fig',n,z(iz));
        
        
        fig = figure;
        set(fig,'name',maintitle,'numbertitle','off')
        
        %First subplot is the magnitude of the whole image
        subplot(1,3,1);
        imagesc(x,y,abs(epsout_ind(:,:,iz,n)));
        title('Magnitude');
        colorbar;
        set(gca,'ydir','normal');

        %Second subplot is the real part of the image
        subplot(1,3,2);
        imagesc(x,y,real(epsout_ind(:,:,iz,n)));       
        title('Real');
        set(gca,'ydir','normal');
        colorbar;

        %third subplot is the imaginary part of the image
        subplot(1,3,3);
        imagesc(x,y,imag(epsout_ind(:,:,iz,n)));
        title('Imaginary');
        colorbar;
        set(gca,'ydir','normal');
        
        set(findall(fig,'-property','FontSize'),'FontSize',15)
        set(fig, 'Position', [100, 100, 1200, 300]);

        % % enable this if you want to save the figures
        if save_ind savefig(filename); end
        if close_ind close(fig); end
        
    end
    
    %same as above: This image being printed is the total combined image
    %using all the s-parameters provided.
    maintitle = sprintf('Combined Image, Z-plane = %.03f %s',z(iz),units);
    filename = sprintf('%s_z%.03f.fig','S_Combo_',z(iz));

    fig = figure;
    set(fig,'name',maintitle,'numbertitle','off')

    subplot(1,3,1);
    imagesc(x,y,abs((epsout_tot(:,:,iz))));
    title('Magnitude');
    colorbar;
    set(gca,'ydir','normal');

    subplot(1,3,2);
    imagesc(x,y,real((epsout_tot(:,:,iz))));
    title('Real');
    set(gca,'ydir','normal');
    colorbar;

    subplot(1,3,3);
    imagesc(x,y,imag((epsout_tot(:,:,iz))));
    title('Imaginary');
    colorbar;
    set(gca,'ydir','normal');
    
    set(findall(fig,'-property','FontSize'),'FontSize',15)
    set(fig, 'Position', [100, 100, 1200, 300]);

    if save_comb savefig(filename); end
    if close_comb close(fig); end
end

toc