% 3D convolution-based SPM employing calibration with dielectric scattering probe
% VERSION DATE: 16/06/2017
% AUTHOR: Denys Shumakov
% in support of Section 4.4.3

% CODE REQUIREMENTS:
% 1. The S-parameter raw data are assumed to be in the (x,y,freq,z) format.
% 2. CO (calibration object) scan area must be of the same size as that of the OUT.

% EXAMPLE
% This code can be run with the provided S-parameter data based on the 
% simulatioin of the F-shape example; see Section 4.3.5.

% DENOISING OPTION
% The CO (calibration object) and the RO (reference object) data can be 
% de-noised for improved image quality. 
% The code below demonstrates the use of de-noising on the CO data.
% To enable de-noising, remove "%" in front of respective line. 
% The de-noising code is provided in the folder "DenoisingCodes".
% For the details of the de-noising algorithm, see [Shumakov, Beaverstone and 
% Nikolova, "De-noising algorithm for enhancing microwave imaging", IET JoE, 2017]
% It is NOT recommended to de-noise the OUT data since this may result in
% loss of information
%%
%% load the data

clear; clc
tic
% load reflection coefficient at port 1 - calibration object (CO) acquisition
load 'S11_CO';
% de-noising CO S11
%[S11_CO_signal,S11_CO_noise]=denoising(S11_CO,0);
S1_CO=S11_CO;
% load reflection coefficient at port 2 - CO acquisition
load 'S22_CO';
% de-noising CO S22
%[S22_CO_signal,S22_CO_noise]=denoising(S22_CO,0);
S2_CO=S22_CO;
% load transmission coefficient from port 2 to port 1 - CO acquisition
load 'S12_CO';
% de-noising CO S12
%[S12_CO_signal,S12_CO_noise]=denoising(S12_CO,0); %enable de-noising algorithm
S3_CO=S12_CO;
% load transmission coefficient from port 1 to port 2 - CO acquisition
load 'S21_CO';
% de-noising CO S21
%[S21_CO_signal,S21_CO_noise]=denoising(S21_CO,0); %enable de-noising algorithm
S4_CO=S21_CO;
% load reflection coefficient at port 1 - object under test (OUT) acquisition
load 'S11_OUT';
S1_OUT=S11_OUT;
% load reflection coefficient at port 2 - OUT acquisition
load 'S22_OUT';
S2_OUT=S22_OUT;
% load transmission coefficient from port 2 to port 1 - OUT acquisition
load 'S12_OUT';
S3_OUT=S12_OUT;
% load transmission coefficient from port 1 to port 2 - OUT acquisition
load 'S21_OUT';
S4_OUT=S21_OUT;
% load reflection coefficient at port 1 - reference object (RO) acquisition
load 'S11_RO';
%[S11_RO_signal,S11_RO_noise]=denoising(S11_RO); %enable de-noising algorithm
S1_RO=S11_RO;
% load reflection coefficient at port 2 - (RO) acquisition
load 'S22_RO';
%[S22_RO_signal,S22_RO_noise]=denoising(S22_RO); %enable de-noising algorithm
S2_RO=S22_RO;
% load transmission coefficient from port 2 to port 1 - RO acquisition
load 'S12_RO';
%[S12_RO_signal,S12_RO_noise]=denoising(S12_RO); %enable de-noising algorithm
S3_RO=S12_RO;
% load transmission coefficient from port 1 to port 2 - RO acquisition
load 'S21_RO';
%[S21_RO_signal,S21_RO_noise]=denoising(S21_RO); %enable de-noising algorithm
S4_RO=S21_RO;


%% initialization

[xnum, ynum, fnum,znum]=size(S1_CO); %xnum, ynum and znum are number of samples along x,y and z; fnum is number of frequency samples
f_start=1; %number of freq sample to start with
f_stop=fnum; %number of freq sample to finish with
f_range=f_start:f_stop; %set the frequency range to be processed
COsize=10; %in mm
scan_step=10; %spatial sampling step in mm
ROcontr=1; %relative permittivity of background (RO)
COcontr=0.1; %known contrast of scattering probe in CO (i.e., the difference between 
% the relative permittivity of the scattering probe and that of the RO)
% Choose type of forward model: Born or Rytov
approx=1; % Born is 1, Rytov is 0. 
% For details, see [Shumakov and Nikolova, "Fast quantitative microwave imaging
% with scattered-power maps", IEEE Trans. MTT, 2017]
scaling_factor=(COsize^2)/(scan_step^2); %calculating scaling factor for the cubic scattering probe. See Eq.(4.102)
co_max=zeros(1,fnum);

%% generating OUT power maps: M_OUT; see Eq.(4.84)

for t=1:znum %loop over all range locations
    for i=f_range %loop over frequency range to be processed
        if approx == 1 %calculating scattered responses within the LQL (or Born) approximation; see Section 4.1
 %subtract incident-field response from tota-field response to get
 %scattered-field response
            deltaS1=S1_OUT(:, :, i)-S1_RO(:, :, i);
            deltaS2=S2_OUT(:, :, i)-S2_RO(:, :, i);
            deltaS3=S3_OUT(:, :, i)-S3_RO(:, :, i);
            deltaS4=S4_OUT(:, :, i)-S4_RO(:, :, i);
            dS1=S1_CO(:, :, i, t)-S1_RO(:, :, i);
            dS2=S2_CO(:, :, i, t)-S2_RO(:, :, i);
            dS3=S3_CO(:, :, i, t)-S3_RO(:, :, i);
            dS4=S2_CO(:, :, i, t)-S4_RO(:, :, i);
        else %calculating scattered responses within the LQLR (or Rytov) approximation
% see [Shumakov and Nikolova, "Fast quantitative microwave imaging
% with scattered-power maps", IEEE Trans. MTT, 2017]
            deltaS1=S1_RO(:, :, i).*log(S1_OUT(:, :, i)./S1_RO(:, :, i));
            deltaS2=S2_RO(:, :, i).*log(S2_OUT(:, :, i)./S2_RO(:, :, i));
            deltaS3=S3_RO(:, :, i).*log(S3_OUT(:, :, i)./S3_RO(:, :, i));
            deltaS4=S4_RO(:, :, i).*log(S4_OUT(:, :, i)./S4_RO(:, :, i));
            dS1=S1_RO(:, :, i).*log(S1_CO(:, :, i, t)./S1_RO(:, :, i));
            dS2=S2_RO(:, :, i).*log(S2_CO(:, :, i, t)./S2_RO(:, :, i));
            dS3=S3_RO(:, :, i).*log(S3_CO(:, :, i, t)./S3_RO(:, :, i));
            dS4=S4_RO(:, :, i).*log(S4_CO(:, :, i, t)./S4_RO(:, :, i));
        end
        
        if t*i==f_start %initialization of arrays for the data in Fourier domain
            S1COF=zeros(xnum,ynum,fnum,znum);
            S2COF=S1COF;
            S3COF=S1COF;
            S4COF=S1COF;
            cmap_out=zeros(xnum,ynum,znum);
            
        end
        
        %2D FFT of the CO and OUT data sets
        S1COF(:,:,i,t)=fft2(dS1);
        S2COF(:,:,i,t)=fft2(dS2);
        S3COF(:,:,i,t)=fft2(dS3);
        S4COF(:,:,i,t)=fft2(dS4);
        S1COUTF=fft2(deltaS1);
        S2COUTF=fft2(deltaS2);
        S3COUTF=fft2(deltaS3);
        S4COUTF=fft2(deltaS4);
        
%calculating the 2D FT of the OUT power maps at each range location and frequency sample
%using all responses
%see Eq.(4.65)
        cmap_out_tmp=S1COUTF.*conj(S1COF(:,:,i,t))+S2COUTF.*conj(S2COF(:,:,i,t))...
            +S3COUTF.*conj(S3COF(:,:,i,t))+S4COUTF.*conj(S4COF(:,:,i,t));
        
        cmap_co_tmp=zeros(xnum,ynum,znum); %initializing the CO power maps array
        for tt=1:znum %loop over all range locations
%calculating the 2D FT of the CO power maps for each range location and each frequency sample
            cmap_co_tmp(:,:,tt)=S1COF(:,:,i,tt).*conj(S1COF(:,:,i,tt))+S2COF(:,:,i,tt).*conj(S2COF(:,:,i,tt))...
                +S3COF(:,:,i,tt).*conj(S3COF(:,:,i,tt))+S4COF(:,:,i,tt).*conj(S4COF(:,:,i,tt));
        end
        
%finding the maximum of the Fourier-domain CO power map over all range locations at each frequency sample
%see Eq.(4.32)
        co_max(:,i)=max(max(max(abs(cmap_co_tmp))));
        cmap_out(:,:,t)=cmap_out(:,:,t)+cmap_out_tmp./co_max(:,i); %frequency combination of the normalized OUT power maps. 
        
    end
end

%% generating CO power maps: M_CO (same as M_OUT but using CO data); see Eq.(4.88)

cmap_co_tmp2=zeros(xnum,ynum,znum); %initializing the CO power maps array
for tt=1:znum %loop over all range locations
    for i=f_range %loop over frequency range to be processed
        
        %calculating the CO power maps for each range location and each frequency sample
        cmap_co_tmp2(:,:,tt)=S1COF(:,:,i,tt).*conj(S1COF(:,:,i,tt))+S2COF(:,:,i,tt).*conj(S2COF(:,:,i,tt))...
            +S3COF(:,:,i,tt).*conj(S3COF(:,:,i,tt))+S4COF(:,:,i,tt).*conj(S4COF(:,:,i,tt));
        
        if tt*i==f_start %initialization of the CO power maps array
            CO_Cmap_Store=zeros(xnum,ynum,znum,znum);
        end
        
        CO_Cmap_Store(:,:,1:end,tt)=CO_Cmap_Store(:,:,1:end,tt)+cmap_co_tmp2(:,:,1:end)./co_max(:,i); %frequency combination of the normalized CO power maps.
        %See Section 4.3.3
    end
end


%% solving linear system of equations at each point in Fourier domain. See Eq.(4.103)

x1=zeros(xnum,ynum,znum); %initializing solution array in Fourier domain
for kx=1:xnum
    for ky=1:ynum
        x1(kx,ky,:)=squeeze(CO_Cmap_Store(kx,ky,:,:))\squeeze(cmap_out(kx,ky,:)); %finding the solution at each kappa in Fourier domain
    end
end

x=zeros(xnum,ynum,znum); %initializing solution array in spatial domain
for t=1:znum
    x(:,:,t)=scaling_factor.*COcontr.*ifftshift(ifft2(x1(:,:,t)))+ROcontr; %finding the solution in spatial domain; see Eq.(4.105)
end

toc

%% Plotting the results

int=0; %interpolation yes 1, no 0. Note that the colormap values may change during certain interpolation settings.
int_rate=0.2; %set interpolation rate

if znum == 1 %plotting 2D images (this section can be also run to get the 3D images slice by slice separately: that is why loop over range locations is still there)
    
    for t=1:znum
        figure('Units', 'centimeters', 'OuterPosition',[0 10 70 14]);
        
        subplot(1, 3, 1) %qualitative outcome
        drawmap=ifftshift(ifft2(cmap_out(:,:,t)));
        
        V=abs(drawmap/max(max(drawmap))); %normalized modulus
        imagesc(0:(xnum-1)*scan_step/10, 0:(ynum-1)*scan_step/10,V);
        colorbar
        colormap parula
        shading interp
        alpha(0.9);
        
        subplot(1, 3, 2) %real part of the relative dielectric permittivity distribution
        drawmap=x(:,:,t);
        if int == 1
            [X,Y]=meshgrid(1:xnum,1:ynum);
            [XI,YI]=meshgrid(1:int_rate:xnum,1:int_rate:ynum);
            V2=interp2(X,Y,drawmap,XI,YI,'spline');
            imagesc(0:(xnum-1)*scan_step/10, 0:(ynum-1)*scan_step/10,real(V2));
        else
            V=real(drawmap);
            imagesc(0:(xnum-1)*scan_step/10, 0:(ynum-1)*scan_step/10,V);
        end
        colorbar
        colormap parula
        shading interp
        alpha(0.9);
        
        subplot(1, 3, 3) %imaginary part of the relative dielectric permittivity distribution
        if int == 1
            [X,Y]=meshgrid(1:xnum,1:ynum);
            [XI,YI]=meshgrid(1:int_rate:xnum,1:int_rate:ynum);
            V2=interp2(X,Y,drawmap,XI,YI,'spline');
            imagesc(0:(xnum-1)*scan_step/10, 0:(ynum-1)*scan_step/10,imag(V2));
        else
            V=imag(drawmap);
            imagesc(0:(xnum-1)*scan_step/10, 0:(ynum-1)*scan_step/10,V);
        end
        colorbar
        colormap parula
        shading interp
        alpha(0.9);
        
        set(gcf, 'color','w');
    end
    
else % plotting 3D images
    
    if int == 1
        Z=ones(xnum,ynum);
        [X,Y]=meshgrid(1:xnum,1:ynum);
        [X1,Y1]=meshgrid(linspace(0,(xnum-1)*scan_step/10,xnum),linspace(0,(ynum-1)*scan_step/10,ynum));
        
        figure
        for k=1:znum
            [XI,YI]=meshgrid(1:int_rate:xnum,1:int_rate:ynum);
            V3(:,:,k)=interp2(X,Y,x(:, :, k), XI,YI);
            surface(X1,Y1,Z.*k,real(V3(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
        
        figure
        for k=1:znum
            [XI,YI]=meshgrid(1:int_rate:xnum,1:int_rate:ynum);
            V3(:,:,k)=interp2(X,Y,x(:, :, k), XI,YI);
            surface(X1,Y1,Z.*k,imag(V3(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
        
    else
        Z=ones(xnum,ynum);
        [X,Y]=meshgrid(1:xnum,1:ynum);
        figure
        for k=1:znum
            surface(X,Y,Z.*k,real(x(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
        
        figure
        for k=1:znum
            surface(X,Y,Z.*k,imag(x(:,:,k)),'FaceColor','texturemap','CDataMapping','scaled')
        end
        colormap parula
        view(3), box on, axis tight square
        shading flat
        set(gca,'YDir','reverse','ZLim',[1 znum])
        colorbar
    end
end