%De-noising algorithm for the complex-valued data as of 16/06/2017 by Denys Shumakov.
%Algorithm works on the real and imaginary parts of the complex data separately.

%CODE REQUIREMENTS:

%1. The S-parameter raw data are assumed to be in the (x,y,freq,z) format.
%2. Param (see line 13) set to 1 - de-noising with convergence criterion enabled. Param set to 0 - de-noising without convergence criterion (maximum number of
%iterations has to be set in line 24).
%3. Interpolation method can be set to "natural" or "nearest" in line 17 (see help for "scatteredInterpolant").
%4. For details of the algorithm, see [Shumakov, Beaverstone and Nikolova "De-noising algorithm for enhancing microwave imaging", IET JoE,
%2017]

function [S_CO_signal,S_CO_noise]=denoising(Sparameter,param)

[xnum, ynum, fnum, znum]=size(Sparameter); %xnum, ynum and znum are number of samples along x,y and z; fnum is number of frequency samples
conn=conndef(2,'minimal'); %setting minimum connectivity region for finding local minima/maxima
method='natural'; %2D interpolation method
Sparameter_signal=zeros(xnum, ynum, fnum,znum);
Sparameter_noise=zeros(xnum, ynum, fnum,znum);
S_CO_signal=zeros(xnum, ynum, fnum,znum);
S_CO_noise=zeros(xnum, ynum, fnum,znum);
conf99=sqrt(2)*erfcinv(2*.01/2); %0.99 probability point of cumulative distribution function of standard normal distribution
upconf1=conf99/sqrt(xnum*ynum); %setting upper confidence bound
max_it=60; %maximum number of iterations

for t=1:znum %loop over all range locations
    for f=1:fnum %loop over all frequency samples
        noi=1; %current number of iterations
        while noi < max_it
            
            if noi == 1 %initialization
                S_CO=Sparameter(:,:,f,t);
            else
                S_CO=Sparameter_noise(:,:,f,t);
            end
            
            reS_CO=squeeze(real(S_CO)); %real part of data
            imS_CO=squeeze(imag(S_CO)); %imag part of data
            
            A=imregionalmin(reS_CO,conn); %finding local minima for real part
            B=imregionalmin(imS_CO,conn); %finding local minima for imag part
            
            [rowa,cola]=find(A); %finding indices of local minima
            [rowb,colb]=find(B); %finding indices of local minima
            
            R=zeros(length(rowa),1);
            P=[rowa cola];
            for j=1:length(rowa)
                R(j)=real(S_CO(rowa(j),cola(j))); %finding values of local minima for real part
            end
            F=scatteredInterpolant(P,R,method); %scattered data interpolation
            
            RI=zeros(length(rowb),1);
            PI=[rowb colb];
            for j=1:length(rowb)
                RI(j)=imag(S_CO(rowb(j),colb(j))); %finding values of local minima for imag part
            end
            FI=scatteredInterpolant(PI,RI,method); %scattered data interpolation
            
            ti = 1:xnum;
            tj=1:ynum;
            [yq,xq] = meshgrid(tj,ti); %creating 2D rectangular grid
            vq_real_min = F(xq,yq); %evaluating interpolating function at meshgrid locations for real part
            vq_imag_min = FI(xq,yq); %evaluating interpolating function at meshgrid locations for imag part
            
            reS_CO=zeros(xnum, ynum);
            imS_CO=zeros(xnum, ynum);
            
            reS_COm=squeeze(real(S_CO)); %real part of data
            imS_COm=squeeze(imag(S_CO)); %imag part of data
            
            AM=imregionalmax(reS_COm,conn); %finding local maxima for real part
            BM=imregionalmax(imS_COm,conn); %finding local maxima for imag part
            
            [rowam,colam]=find(AM); %finding indices of local maxima for real part
            [rowbm,colbm]=find(BM); %finding indices of local maxima for imag part
            
            RM=zeros(length(rowam),1);
            PM=[rowam colam];
            for j=1:length(rowam)
                RM(j)=real(S_CO(rowam(j),colam(j))); %finding values of local maxima for real part
            end
            FM=scatteredInterpolant(PM,RM,method);
            
            RIM=zeros(length(rowbm),1);
            PIM=[rowbm colbm];
            for j=1:length(rowbm)
                RIM(j)=imag(S_CO(rowbm(j),colbm(j))); %finding values of local maxima for imag part
            end
            FIM=scatteredInterpolant(PIM,RIM,method); %scattered data interpolation
            
            ti = 1:xnum;
            tj=1:ynum;
            [yq,xq] = meshgrid(tj,ti); %creating 2D rectangular grid
            vq_real_max = FM(xq,yq); %evaluating interpolating function at meshgrid locations for real part
            vq_imag_max = FIM(xq,yq); %evaluating interpolating function at meshgrid locations for imag part
            
            reS_COm=zeros(xnum, ynum);
            imS_COm=zeros(xnum, ynum);
            
            
            mean_real=(vq_real_max+vq_real_min)/2; %finding the arithmetic mean of two surfaces for real part
            mean_imag=(vq_imag_max+vq_imag_min)/2; %finding the arithmetic mean of two surfaces for imag part
            trend=mean_real+1i*mean_imag; %finding the trend
            
            reS_CO_new=real(S_CO)-mean_real;
            imS_CO_new=imag(S_CO)-mean_imag;
            residual=reS_CO_new+1i*imS_CO_new; %finding the residual
            
            
            Sparameter_signal(:,:,f,t)=Sparameter_signal(:,:,f,t)+trend; %signal component
            Sparameter_noise(:,:,f,t)=residual; %resultant noise component
            
            if param == 1 %convergence check
                signal=reshape(Sparameter_signal(:,:,f,t),[],1);
                noise2=reshape(Sparameter_noise(:,:,f,t),[],1);
                
                var1=max(abs(xcorr(signal,noise2,'coeff'))); %cross-correlation between the signal and noise
                
                if var1 <= upconf1                    
                    break                    
                else
                    if noi == max_it-1
                        disp('convergence criterion not reached for f_sample ='); %notification in case algorithm did not converge
                        disp(f)
                    end
                    noi=noi+1;
                end
            else
                noi=noi+1;
            end
        end
        
        S_CO_signal(:,:,f,t)=Sparameter_signal(:,:,f,t); %storing the signal at each frequency and range location
        S_CO_noise(:,:,f,t)=Sparameter_noise(:,:,f,t); %storing the noise at each frequency and range location
    end
end

