calculate_cs_uncertainty.m

function [Ux, Uy] = calculate_cs_uncertainty(im1, im2, x1, y1, CovDCx, CovDCy, window_size, window_resolution)
%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

    % 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, :);

    %% 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);
    L = size(CovDCx);

    % calculate product of shifted images to evaluate Correlation function C
    C0_full = im1 .* im2;
    Cp1x_full = im1 .* im2shiftx;
    Cn1x_full = im1shiftx .* im2;
    Cp1y_full = im1 .* im2shifty;
    Cn1y_full = im1shifty .* im2;

    %Summing the correlation functions
    C0 = sum(C0_full(:));
    Cp1x = sum(Cp1x_full(:));
    Cn1x = sum(Cn1x_full(:));
    Cp1y = sum(Cp1y_full(:));
    Cn1y = sum(Cn1y_full(:));

    % Coorelation function positive and negative x and y values as mentione din
    % the paper
    Cpnx = (Cp1x + Cn1x)/2;
    Cpny = (Cp1y + Cn1y)/2;

    % calculate extents of present image
    xmin1 = x1 - ceil(Nx/2)+1;
    xmax1 = x1 + floor(Nx/2);
    ymin1 = y1 - ceil(Ny/2)+1;
    ymax1 = y1 + floor(Ny/2);

    % extract covariance 
    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]));
    
    %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 = (1/2) * ( log(Cpnx + SigDCx/2) - log(Cpnx - SigDCx/2) ) / (2*log(C0) - log(Cpnx + SigDCx/2) - log(Cpnx - SigDCx/2));
    Uy = (1/2) * ( log(Cpny + SigDCy/2) - log(Cpny - SigDCy/2) ) / (2*log(C0) - log(Cpny + SigDCy/2) - log(Cpny - SigDCy/2));

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