correlation_statistics.m

function [Ux,Uy]=correlation_statistics(im1, im2, window_size, window_resolution, X, Y, dx, dy)
%This function calculates the uncertainty in X and Y displacement by calculating
% the variance DeltaC (difference in cross correlation function)given a pair of 
% corrected images(with estimated displacement field in Iterative window_size Deformation) 

% This is the Correlation Statistics Method developed by B. Wieneke (2014)

%Input
% im1 and im2 are deformed interrogation images after last pass
% window_size: window size for processing
% window_resolution: window resolution for processing
% X,Y the grid points for velocity uncertainty evaluation
% dx,dy: diameter region over which autocovariance sum is to be evaluated
% default. dx = 2, dy = 2
%Output
% Ux,Uy: Uncertainty in X and Y direction for each grid point

% The equation to be used for uncertainty in displacement
% Sigu=(dx/2)*(log (Cpn + SigDC/2) - log (Cpn - SigDC/2))/(2logC0 - log (Cpn + SigDC/2) - log (Cpn - SigDC/2));

% So we have to calculate SigDCx and SigDCy which are essentially variance
% of DC or DeltaC which is difference in correlation function defined by equation 
% 4 in the paper.

%In terms of implmentation two paths can be followed. Evaluation the
%covariance function for each window or for the whole image and then
%finding the sum and corresponding uncertainty for each window. The latter
%is computationally efficient. The former has been tried before but this
%code is written for the whole image covariance calculation

%The evaluation requires a guassian recursive filter which has been
%implemented following the cited paper. However, the smoothing and
%filtering required in this algorithm may need to be tuned in future for
%optimal performance and I have seen some variation in results for
%different values of this parameters. The basic structure of the algorithm
%is fine.

%This is implemented by Sayantan Bhattacharya
% Modified by Lalit Rajendran, Nov. 2019

    % set default values for summation neighborhood if not specified
    if nargin < 7
        dx = 2; %4;
        dy = 2; %4;
    end
    %% Find SigDCx and SigDCy

    % mean subtraction on images
    im1 = im1 - mean(im1(:));
    im2 = im2 - mean(im2(:));

    % Shift images in x and y direction by 1 pixel.
    im1shiftx=zeros(size(im1));
    im1shiftx(:,1:end-1)=im1(:,2:end);

    im2shiftx=zeros(size(im2));
    im2shiftx(:,1:end-1)=im2(:,2:end);

    im1shifty=zeros(size(im1));
    im1shifty(1:end-1,:)=im1(2:end,:);

    im2shifty=zeros(size(im2));
    im2shifty(1:end-1,:)=im2(2:end,:);

    % calculate product of shifted images to evaluate Correlation function C
    C0reg=im1.*im2;
    Cp1xreg=im1.*im2shiftx;
    Cn1xreg=im1shiftx.*im2;
    Cp1yreg=im1.*im2shifty;
    Cn1yreg=im1shifty.*im2;

    %Calculate delta of correlation function or Delta C or DC
    %Calculate elemental DC matrix DCix and DCiy given by equation 7 in x and y
    DCix=im1.*im2shiftx - im1shiftx.*im2;
    DCiy=im1.*im2shifty - im1shifty.*im2;

    [Sy,Sx]=size(DCix);

    % Apply 1-2-1/4 smoothing filter on DCix and DCiy in x and y direction
    % respectively to get proper covariance matrix
    DCix=filter([1 2 1],4,DCix,[],2);
    DCiy=filter([1 2 1],4,DCiy,[],1);

%     % mean subtraction (lalit)
%     DCix = DCix - mean(DCix(:), 'omitnan');
%     DCiy = DCiy - mean(DCiy(:), 'omitnan');
    
    % Pad the array with 4 rows and columns of zeros for calculating
    % autocovariance matrix
    DCix=padarray(DCix,[4 4],0);
    DCiy=padarray(DCiy,[4 4],0);

    %Initialize covariance matrix
    covDcx=zeros(size(DCix));
    covDcy=zeros(size(DCiy));

    count=1;
    %Loop over dx-by-dy subregions to calculate covariance in DeltaC function
    for j=1+dx:Sx+dx
        for i=1+dy:Sy+dy

            % Select subregion of over +- 4 pixels to evaluate cross covariance
            % terms
            subDcix=DCix(i-dy:i+dy,j-dx:j+dx);
            subDciy=DCiy(i-dy:i+dy,j-dx:j+dx);

            %Determining Autocovariance using xcorr2
            Sautox=xcorr2(subDcix);
            Sautoy=xcorr2(subDciy);

            Lx=size(Sautox,1);
            px=(Lx+1)/2;
            S00x=Sautox(px,px);
            Sautox=Sautox(px-dy:px+dy,px-dx:px+dx);
            qx=(px+1)/2;
            %normalizing the autocovariance matrix
            Snormx=Sautox./S00x;
            % Finding normalized values within 5%
            yind1=min(abs(find(Snormx(:,qx)<0.05)-qx));
            xind1=min(abs(find(Snormx(qx,:)<0.05)-qx));
            % Summing the covariance matrix for only those values
            Sumautox=sum(sum(Sautox(qx-yind1+1:qx+yind1-1,qx-xind1+1:qx+xind1-1)));
            covDcx(i,j)=Sumautox;

            Ly=size(Sautoy,1);
            py=(Ly+1)/2;
            S00y=Sautoy(py,py);
            Sautoy=Sautoy(py-dy:py+dy,py-dx:py+dx);
            qy=(py+1)/2;
            %normalizing the autocovariance matrix
            Snormy=Sautoy./S00y;
            % Finding normalized values within 5%
            yind2=min(abs(find(Snormy(:,qy)<0.05)-qy));
            xind2=min(abs(find(Snormy(qy,:)<0.05)-qy));
            % Summing the covariance matrix for only those values
            Sumautoy=sum(sum(Sautoy(qy-yind2+1:qy+yind2-1,qy-xind2+1:qy+xind2-1)));
            covDcy(i,j)=Sumautoy;

            count=count+1;
    %         keyboard;
        end
    end

    % Eliminationg the border values
    covDcx=covDcx(1+dy:Sy+dy,1+dx:Sx+dx);
    covDcy=covDcy(1+dy:Sy+dy,1+dx:Sx+dx);

    CovDCx = covDcx;
    CovDCy = covDcy;
    %% finding window_sizes from the whole image and then evaluating uncertainty propagation
    % for each grid point

    Nx = window_size(1);
    Ny = window_size(2);
    L=size(im1);
    X=X(:);
    Y=Y(:);
    Ux=zeros(size(X));
    Uy=zeros(size(Y));

    % window_size masking filter
    sfilt1 = windowmask([Nx Ny],[window_resolution(1, 1) window_resolution(1, 2)]);
    sfilt2 = windowmask([Nx Ny],[window_resolution(2, 1) window_resolution(2, 2)]);

    for n=1:length(X)
        x1 = X(n);
        y1 = Y(n);
        % x2 = X(n);
        % y2 = Y(n);

        xmin1 = x1- ceil(Nx/2)+1;
        xmax1 = x1+floor(Nx/2);
        % xmin2 = x2- ceil(Nx/2)+1;
        % xmax2 = x2+floor(Nx/2);
        ymin1 = y1- ceil(Ny/2)+1;
        ymax1 = y1+floor(Ny/2);
        % ymin2 = y2- ceil(Ny/2)+1;
        % ymax2 = y2+floor(Ny/2);


        %Given grid points and window sizes evaluate the correlation functions 
        %for the particluar window sizes
        if xmin1<1 || ymin1<1 || xmax1>Nx || ymax1>Ny
            %For boundary points
            C0t( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) )   = C0reg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            Cp1xt( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = Cp1xreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            Cn1xt( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = Cn1xreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));

            Cp1yt( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = Cp1yreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            Cn1yt( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = Cn1yreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));

            CovDCxt( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) )=CovDCx( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            CovDCyt( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) )=CovDCy( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
        else
            C0t =C0reg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            Cp1xt=Cp1xreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            Cn1xt=Cn1xreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));

            Cp1yt =Cp1yreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            Cn1yt= Cn1yreg( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));

            CovDCxt=CovDCx( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
            CovDCyt=CovDCy( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));


        end

        % If spatial filtering is used these lines should be uncommented
        C0t=C0t.*sfilt1;
        Cp1xt=Cp1xt.*sfilt1;
        Cn1xt=Cn1xt.*sfilt1;
        Cp1yt=Cp1yt.*sfilt2;
        Cn1yt=Cn1yt.*sfilt2;
        CovDCxt=CovDCxt.*sfilt1;
        CovDCyt=CovDCyt.*sfilt2;

        %Summing the correlation functions
        C0=sum(C0t(:));
        Cp1x=sum(Cp1xt(:));
        Cn1x=sum(Cn1xt(:));
        Cp1y=sum(Cp1yt(:));
        Cn1y=sum(Cn1yt(:));
        
        % Coorelation function positive and negative x and y values as mentione din
        % the paper
        Cpnx=(Cp1x+Cn1x)/2;
        Cpny=(Cp1y+Cn1y)/2;

        %Sigma of Delta Correlation function in x and y direction
        SigDCx=sqrt(sum(CovDCxt(:)));
        SigDCy=sqrt(sum(CovDCyt(:)));

        %% Calculate uncertainties Ux and Uy
        %Using uncertainty propagation equation
        Ux(n)=(1/2)*(log (Cpnx + SigDCx/2) - log (Cpnx - SigDCx/2))/(2*log(C0) - log (Cpnx + SigDCx/2) -log (Cpnx - SigDCx/2));
        Uy(n)=(1/2)*(log (Cpny + SigDCy/2) - log (Cpny - SigDCy/2))/(2*log(C0) - log (Cpny + SigDCy/2) -log (Cpny - SigDCy/2));
    end

    % only retain real part
    Ux = real(Ux);
    Uy = real(Uy);
end

function [Gout]=gauss_recursive_filter(Gin,winres,dim)
    % This is a fast gaussian recursive filter following Lukin 2007, Young 1995

    % Input window_size resolution (winres), 1d or 2d filter implementation (dim=1
    % or 2) and Gin matrix or vector that needs to be filtered and output will
    % be Gout

    filtrad=winres/8;
    if filtrad>=2.5
        qFactor = 0.98711*filtrad - 0.96330;          
    elseif filtrad>=0.5 && filtrad<2.5
        qFactor = 3.97156 -4.14554*sqrt(1-0.26891*filtrad);
    end
    b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
    b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
    b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
    b3Coeff = 0.422205 * qFactor * qFactor * qFactor;
    normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);
    vDenCoeff = [b0Coeff, -b1Coeff, -b2Coeff, -b3Coeff] / b0Coeff;

    if dim==2

        tempg = filter(normalizationCoeff, vDenCoeff, Gin,[],2);
        tempg1 = filter(normalizationCoeff, vDenCoeff, tempg(:,end:-1:1),[],2);
        tempg2 = filter(normalizationCoeff, vDenCoeff, tempg1,[],1);
        Gout = filter(normalizationCoeff, vDenCoeff, tempg2(end:-1:1,:),[],1);
        Gout=(filtrad*sqrt(2*pi))^2.*Gout(end:-1:1,end:-1:1);

    elseif dim==1
        tempg = filter(normalizationCoeff, vDenCoeff, Gin);
        Gout= filter(normalizationCoeff, vDenCoeff, tempg(end:-1:1));
        Gout=(filtrad*sqrt(2*pi)).*Gout(end:-1:1);
    end

end