calculate_covariance_matrix.m

function [CovDCx, CovDCy] = calculate_covariance_matrix(im1, im2, 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 < 3
        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 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;
        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;

end