Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

prana/PIVwindowed.m

Go to file
 
 
Cannot retrieve contributors at this time
executable file 835 lines (715 sloc) 34.8 KB
Raw Blame
Open in GitHub Desktop
function [X,Y,U,V,C,Dia,Corrplanes,uncertainty2D,SNRmetric]=PIVwindowed(im1,im2,tcorr,window,res,zpad,D,Zeromean,Peaklocator,Peakswitch,fracval,saveplane,X,Y,uncertainty,e,Uin,Vin)
% --- DPIV Correlation ---
imClass = 'double';
% This file is part of prana, an open-source GUI-driven program for
% calculating velocity fields using PIV or PTV.
%
% Copyright (C) 2012 Virginia Polytechnic Institute and State
% University
%
% Copyright 2014. Los Alamos National Security, LLC. This material was
% produced under U.S. Government contract DE-AC52-06NA25396 for Los
% Alamos National Laboratory (LANL), which is operated by Los Alamos
% National Security, LLC for the U.S. Department of Energy. The U.S.
% Government has rights to use, reproduce, and distribute this software.
% NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY
% WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF
% THIS SOFTWARE. If software is modified to produce derivative works,
% such modified software should be clearly marked, so as not to confuse
% it with the version available from LANL.
%
% prana is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%convert input parameters
im1=cast(im1,imClass);
im2=cast(im2,imClass);
L=size(im1);
%convert to gridpoint list
X=X(:);
Y=Y(:);
% keyboard;
%preallocate velocity fields and grid format
Nx = window(1);
Ny = window(2);
if nargin <=17
Uin = zeros(length(X),1,imClass);
Vin = zeros(length(X),1,imClass);
end
if Peakswitch
Uin=repmat(Uin(:,1),[1 3]);
Vin=repmat(Vin(:,1),[1 3]);
U = zeros(length(X),3,imClass);
V = zeros(length(X),3,imClass);
C = zeros(length(X),3,imClass);
Dia = zeros(length(X),3,imClass);
else
U = zeros(length(X),1,imClass);
V = zeros(length(X),1,imClass);
C = zeros(length(X),1,imClass);
Dia = zeros(length(X),1,imClass);
end
%sets up extended domain size
if zpad~=0
Sy=2*Ny;
Sx=2*Nx;
elseif strcmpi(tcorr,'DCC')
Sy = res(1,2)+res(2,2)-1;
Sx = res(1,1)+res(2,1)-1;
else
Sy=Ny;
Sx=Nx;
end
%fftshift indicies
fftindy = [ceil(Sy/2)+1:Sy 1:ceil(Sy/2)];
fftindx = [ceil(Sx/2)+1:Sx 1:ceil(Sx/2)];
%window masking filter
sfilt1 = windowmask([Sx Sy],[res(1, 1) res(1, 2)]);
sfilt2 = windowmask([Sx Sy],[res(2, 1) res(2, 2)]);
% sfilt12 = ifft2(fft2(sfilt2).*conj(fft2(sfilt1)));
%
% keyboard
%correlation plane normalization function (always off)
cnorm = ones(Ny,Nx,imClass);
% s1 = fftn(sfilt1,[Sy Sx]);
% s2 = fftn(sfilt2,[Sy Sx]);
% S21 = s2.*conj(s1);
%
% %Standard Fourier Based Cross-Correlation
% iS21 = ifftn(S21,'symmetric');
% iS21 = iS21(fftindy,fftindx);
% cnorm = 1./iS21;
% cnorm(isinf(cnorm)) = 0;
%RPC spectral energy filter
spectral = fftshift(energyfilt(Sx,Sy,D,0));
% This is a check for the fractionally weighted correlation. We won't use
% the spectral filter with FWC or GCC
if strcmpi(tcorr,'FWC')
frac = fracval;
spectral = ones(size(spectral));
elseif strcmpi(tcorr,'GCC')
frac = 1;
spectral = ones(size(spectral));
else
frac = 1;
end
% For dynamic rpc flip this switch which allows for dynamic calcuation of
% the spectral function using the diameter of the autocorrelation.
if strcmpi(tcorr,'DRPC')
dyn_rpc = 1;
else
dyn_rpc = 0;
end
if saveplane
Corrplanes=zeros(Sy,Sx,length(X),imClass);
else
Corrplanes = 0;
end
% Intializing SBRmetric and uncertainty 2D structure
SNRmetric.PPR=zeros(length(X),1);
SNRmetric.MI=zeros(length(X),1);
uncertainty2D.Upprx=zeros(length(X),1);
uncertainty2D.Uppry=zeros(length(X),1);
uncertainty2D.UpprxLB=zeros(length(X),1);
uncertainty2D.UppryLB=zeros(length(X),1);
uncertainty2D.UpprxUB=zeros(length(X),1);
uncertainty2D.UppryUB=zeros(length(X),1);
uncertainty2D.UmixLB=zeros(length(X),1);
uncertainty2D.UmiyLB=zeros(length(X),1);
uncertainty2D.UmixUB=zeros(length(X),1);
uncertainty2D.UmiyUB=zeros(length(X),1);
uncertainty2D.Autod=zeros(length(X),1);
uncertainty2D.Ixx=zeros(length(X),1);
uncertainty2D.Iyy=zeros(length(X),1);
uncertainty2D.biasx=zeros(length(X),1);
uncertainty2D.biasy=zeros(length(X),1);
uncertainty2D.Neff=zeros(length(X),1);
uncertainty2D.Uimx=zeros(length(X),1);
uncertainty2D.Uimy=zeros(length(X),1);
uncertainty2D.Nump=zeros(length(X),1);
uncertainty2D.Ucsx=zeros(length(X),1);
uncertainty2D.Ucsy=zeros(length(X),1);
switch upper(tcorr)
%Standard Cross Correlation
case 'SCC'
if size(im1,3) == 3
Gens=zeros(Sy,Sx,3,imClass);
for n=1:length(X)
%apply the second order discrete window offset
x1 = X(n) - floor(round(Uin(n))/2);
x2 = X(n) + ceil(round(Uin(n))/2);
y1 = Y(n) - floor(round(Vin(n))/2);
y2 = Y(n) + ceil(round(Vin(n))/2);
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);
for r=1:size(im1,3)
%find the image windows
zone1 = im1( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]),r);
zone2 = im2( max([1 ymin2]):min([L(1) ymax2]),max([1 xmin2]):min([L(2) xmax2]),r);
if size(zone1,1)~=Ny || size(zone1,2)~=Nx
w1 = zeros(Ny,Nx);
w1( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = zone1;
zone1 = w1;
end
if size(zone2,1)~=Ny || size(zone2,2)~=Nx
w2 = zeros(Ny,Nx);
w2( 1+max([0 1-ymin2]):Ny-max([0 ymax2-L(1)]),1+max([0 1-xmin2]):Nx-max([0 xmax2-L(2)]) ) = zone2;
zone2 = w2;
end
if Zeromean==1
zone1=zone1-mean(mean(zone1));
zone2=zone2-mean(mean(zone2));
end
%apply the image spatial filter
region1 = (zone1).*sfilt1;
region2 = (zone2).*sfilt2;
%FFTs and Cross-Correlation
f1 = fftn(region1,[Sy Sx]);
f2 = fftn(region2,[Sy Sx]);
P21 = f2.*conj(f1);
%Standard Fourier Based Cross-Correlation
G = ifftn(P21,'symmetric');
G = G(fftindy,fftindx);
Gens(:,:,r) = abs(G);
end
G = mean(Gens,3);
%subpixel estimation
[U(n,:),V(n,:),Ctemp,Dtemp]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D);
if Peakswitch
C(n,:)=Ctemp;
Dia(n,:)=Dtemp;
end
if saveplane
Corrplanes(:,:,n) = G;
end
end
else
for n=1:length(X)
%apply the second order discrete window offset
x1 = X(n) - floor(round(Uin(n))/2);
x2 = X(n) + ceil(round(Uin(n))/2);
y1 = Y(n) - floor(round(Vin(n))/2);
y2 = Y(n) + ceil(round(Vin(n))/2);
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);
%find the image windows
zone1 = im1( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
zone2 = im2( max([1 ymin2]):min([L(1) ymax2]),max([1 xmin2]):min([L(2) xmax2]));
if size(zone1,1)~=Ny || size(zone1,2)~=Nx
w1 = zeros(Ny,Nx);
w1( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = zone1;
zone1 = w1;
end
if size(zone2,1)~=Ny || size(zone2,2)~=Nx
w2 = zeros(Ny,Nx);
w2( 1+max([0 1-ymin2]):Ny-max([0 ymax2-L(1)]),1+max([0 1-xmin2]):Nx-max([0 xmax2-L(2)]) ) = zone2;
zone2 = w2;
end
if Zeromean==1
zone1=zone1-mean(mean(zone1));
zone2=zone2-mean(mean(zone2));
end
%apply the image spatial filter
region1 = (zone1).*sfilt1;
region2 = (zone2).*sfilt2;
%FFTs and Cross-Correlation
f1 = fftn(region1,[Sy Sx]);
f2 = fftn(region2,[Sy Sx]);
P21 = f2.*conj(f1);
%Standard Fourier Based Cross-Correlation
G = ifftn(P21,'symmetric');
G = G(fftindy,fftindx);
G = abs(G);
% Minimum value of the correlation plane is subtracted
G=G-min(G(:));
%subpixel estimation
[U(n,:),V(n,:),Ctemp,Dtemp,DXtemp,DYtemp]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D);
if Peakswitch
C(n,:)=Ctemp;
Dia(n,:)=Dtemp;
end
if saveplane
Corrplanes(:,:,n) = G;
end
% Evaluate uncertainty options for SCC
if uncertainty.ppruncertainty(e)==1
%SNR calculation the other output arguments of Cal_SNR
%are Maximum peak value,PRMSR,PCE,ENTROPY
metric='PPR';
PPRval = Cal_SNR(G,metric);
% Save the SNR metrics
SNRmetric.PPR(n)=PPRval;
% Evaluate PPR Uncertainty
% John J Charonko Model
[Ux,Uy,~,~,~,~]=calibration_based_uncertainty('PPR_Charonkomodel',PPRval,upper(tcorr));
uncertainty2D.Upprx(n)=Ux;
uncertainty2D.Uppry(n)=Uy;
% Xue Model
[~,~,UxLB,UxUB,UyLB,UyUB]=calibration_based_uncertainty('PPR_Xuemodel',PPRval,upper(tcorr));
uncertainty2D.UpprxLB(n)=UxLB;
uncertainty2D.UppryLB(n)=UyLB;
uncertainty2D.UpprxUB(n)=UxUB;
uncertainty2D.UppryUB(n)=UyUB;
end
if uncertainty.miuncertainty(e)==1
%Autocorrelations
P11 = f1.*conj(f1);
P22 = f2.*conj(f2);
Auto1 = ifftn(P11,'symmetric');
Auto2 = ifftn(P22,'symmetric');
Auto1 = Auto1(fftindy,fftindx);
Auto2 = Auto2(fftindy,fftindx);
Auto1=abs(Auto1);
Auto2=abs(Auto2);
nAuto1 = Auto1-min(Auto1(:)); % Autocorrelation plane of image 1
nAuto2 = Auto2-min(Auto2(:)); % Autocorrelation plane of image 2
% 3 pt Gaussian fit to Autocorrelation Diameter
[~,~,~,~,Dauto1x3,Dauto1y3,~]=subpixel(nAuto1,Sx,Sy,cnorm,1,0,D);
[~,~,~,~,Dauto2x3,Dauto2y3,~]=subpixel(nAuto2,Sx,Sy,cnorm,1,0,D);
Diap1=sqrt(Dauto1x3*Dauto1y3/2);
Diap2=sqrt(Dauto2x3*Dauto2y3/2);
% Average Autocorrelation Diameter
Autod=mean([Diap1 Diap2]);
uncertainty2D.Autod(n)=Autod;
%MI Calculation
INTS1 = max(region1(:));
INTS2 = max(region2(:));
[MIval,~,~,~,~,~] = MI_Cal_SCC(G,nAuto1,nAuto2,INTS1,INTS2,Dauto1x3,Dauto1y3,Dauto2x3,Dauto2y3,Sx,Sy,fftindx,fftindy);
SNRmetric.MI(n)=MIval;
% Estimate MI uncertainty
[~,~,UxLB,UxUB,UyLB,UyUB]=calibration_based_uncertainty('MI_Xuemodel',MIval,upper(tcorr));
uncertainty2D.UmixLB(n)=UxLB;
uncertainty2D.UmiyLB(n)=UyLB;
uncertainty2D.UmixUB(n)=UxUB;
uncertainty2D.UmiyUB(n)=UyUB;
end
if uncertainty.mcuncertainty(e)==1
if uncertainty.miuncertainty(e)==1
MIest=SNRmetric.MI(n);
[Ixx,Iyy,biasx,biasy,Neff,~]=Moment_of_correlation(P21,f1,f2,Sx,Sy,cnorm,D,fftindx,fftindy,G,DXtemp,DYtemp,region1,region2,MIest);
else
MIest=-1;
[Ixx,Iyy,biasx,biasy,Neff,Autod]=Moment_of_correlation(P21,f1,f2,Sx,Sy,cnorm,D,fftindx,fftindy,G,DXtemp,DYtemp,region1,region2,MIest);
uncertainty2D.Autod(n)=Autod;
end
uncertainty2D.Ixx(n)=Ixx;
uncertainty2D.Iyy(n)=Iyy;
uncertainty2D.biasx(n)=biasx;
uncertainty2D.biasy(n)=biasy;
uncertainty2D.Neff(n)=Neff;
end
% % perform image matching uncertainty calculation (moved here by lalit)
% if uncertainty.imuncertainty(e)==1
% % Perform image matching uncertainty estimation
% [deltax,deltay,Npartu,Npartv] = original_particle_disparity(region1, region2, Sx/2, Sy/2, unique(res));
% % keyboard;
% % assign values to the uncertainty structure
% uncertainty2D.Uimx(n) = deltax;
% uncertainty2D.Uimy(n) = deltay;
% Npart(:,:,1) = Npartu;
% Npart(:,:,2) = Npartv;
% uncertainty2D.Nump(n) = mean(Npart,3);
% clear Npart
% end
% %% calculate uncertainty using correlation statistics (lalit)
% if uncertainty.csuncertainty(e) == 1
% [sigma_cs_x, sigma_cs_y] = correlation_statistics(region1, region2, window, res, Sx/2, Sy/2);
% % [sigma_cs_x, sigma_cs_y] = correlation_statistics(im1d, im2d, Wsize(e,:),Wres(:, :, e), X(Eval>=0),Y(Eval>=0));
% % [sigma_cs_x, sigma_cs_y] = correlation_statistics(im1, im2, Wsize(e,:),Wres(:, :, e), X(Eval>=0),Y(Eval>=0));
% % keyboard;
% % copy cs uncertainties in masked regions
% uncertainty2D.Ucsx(n) = sigma_cs_x;
% uncertainty2D.Ucsy(n) = sigma_cs_y;
% end
end
end
%Direct Cross Correlation
case 'DCC'
%initialize correlation tensor
CC = zeros(Sy,Sx,length(X),imClass);
if size(im1,3) == 3
Gens=zeros(res(1,2)+res(2,2)-1,res(1,1)+res(2,1)-1,3,imClass);
for n=1:length(X)
%apply the second order discrete window offset
x1 = X(n) - floor(round(Uin(n))/2);
x2 = X(n) + ceil(round(Uin(n))/2);
y1 = Y(n) - floor(round(Vin(n))/2);
y2 = Y(n) + ceil(round(Vin(n))/2);
xmin1 = x1- ceil(res(1,1)/2)+1;
xmax1 = x1+floor(res(1,1)/2);
xmin2 = x2- ceil(res(2,1)/2)+1;
xmax2 = x2+floor(res(2,1)/2);
ymin1 = y1- ceil(res(1,2)/2)+1;
ymax1 = y1+floor(res(1,2)/2);
ymin2 = y2- ceil(res(2,2)/2)+1;
ymax2 = y2+floor(res(2,2)/2);
for r=1:size(im1,3)
%find the image windows
zone1 = im1( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]),r );
zone2 = im2( max([1 ymin2]):min([L(1) ymax2]),max([1 xmin2]):min([L(2) xmax2]),r );
if size(zone1,1)~=Ny || size(zone1,2)~=Nx
w1 = zeros(res(1,2),res(1,1));
w1( 1+max([0 1-ymin1]):res(1,2)-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):res(1,1)-max([0 xmax1-L(2)]) ) = zone1;
zone1 = w1;
end
if size(zone2,1)~=Ny || size(zone2,2)~=Nx
w2 = zeros(res(2,2),res(2,1));
w2( 1+max([0 1-ymin2]):res(2,2)-max([0 ymax2-L(1)]),1+max([0 1-xmin2]):res(2,1)-max([0 xmax2-L(2)]) ) = zone2;
zone2 = w2;
end
if Zeromean==1
zone1=zone1-mean(mean(zone1));
zone2=zone2-mean(mean(zone2));
end
%apply the image spatial filter
region1 = (zone1);%.*sfilt1;
region2 = (zone2);%.*sfilt2;
%correlation done using xcorr2 which is faster then fft's
%for strongly uneven windows.
%G = xcorr2(region2,region1);
% This is a stripped out version of xcorr2
G = conv2(region2, rot90(conj(region1),2));
region1_std = std(region1(:));
region2_std = std(region2(:));
if region1_std == 0 || region2_std == 0
G = zeros(Sy,Sx);
else
G = G/std(region1(:))/std(region2(:))/length(region1(:));
end
Gens(:,:,r) = G;
%store correlation matrix
end
G = mean(Gens,3);
%subpixel estimation
[U(n,:),V(n,:),Ctemp,Dtemp]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D);
if Peakswitch
C(n,:)=Ctemp;
Dia(n,:)=Dtemp;
end
if saveplane
Corrplanes(:,:,n) = G;
end
end
else
for n=1:length(X)
%apply the second order discrete window offset
x1 = X(n) - floor(round(Uin(n))/2);
x2 = X(n) + ceil(round(Uin(n))/2);
y1 = Y(n) - floor(round(Vin(n))/2);
y2 = Y(n) + ceil(round(Vin(n))/2);
xmin1 = x1- ceil(res(1,1)/2)+1;
xmax1 = x1+floor(res(1,1)/2);
xmin2 = x2- ceil(res(2,1)/2)+1;
xmax2 = x2+floor(res(2,1)/2);
ymin1 = y1- ceil(res(1,2)/2)+1;
ymax1 = y1+floor(res(1,2)/2);
ymin2 = y2- ceil(res(2,2)/2)+1;
ymax2 = y2+floor(res(2,2)/2);
%find the image windows
zone1 = im1( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]) );
zone2 = im2( max([1 ymin2]):min([L(1) ymax2]),max([1 xmin2]):min([L(2) xmax2]) );
if size(zone1,1)~=res(1,2) || size(zone1,2)~=res(1,1)
w1 = zeros(res(1,2),res(1,1));
w1( 1+max([0 1-ymin1]):res(1,2)-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):res(1,1)-max([0 xmax1-L(2)]) ) = zone1;
zone1 = w1;
end
if size(zone2,1)~=res(2,2) || size(zone2,2)~=res(2,1)
w2 = zeros(res(2,2),res(2,1));
w2( 1+max([0 1-ymin2]):res(2,2)-max([0 ymax2-L(1)]),1+max([0 1-xmin2]):res(2,1)-max([0 xmax2-L(2)]) ) = zone2;
zone2 = w2;
end
if Zeromean==1
zone1=zone1-mean(zone1(:));
zone2=zone2-mean(zone2(:));
end
%apply the image spatial filter
region1 = (zone1);%.*sfilt1;
region2 = (zone2);%.*sfilt2;
%correlation done using xcorr2 which is faster then fft's
%for strongly uneven windows.
%G = xcorr2(region2,region1);
% This is a stripped out version of xcorr2
G = conv2(region2, rot90(conj(region1),2));
region1_std = std(region1(:));
region2_std = std(region2(:));
if region1_std == 0 || region2_std == 0
G = zeros(Sy,Sx);
else
G = G/std(region1(:))/std(region2(:))/length(region1(:));
end
%subpixel estimation
[U(n,:),V(n,:),Ctemp,Dtemp]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D);
if Peakswitch
C(n,:)=Ctemp;
Dia(n,:)=Dtemp;
end
if saveplane
Corrplanes(:,:,n) = G;
end
end
end
%Robust Phase Correlation
case {'RPC','DRPC','GCC','FWC'}
if size(im1,3) == 3
Gens=zeros(Sy,Sx,3,imClass);
for n=1:length(X)
%apply the second order discrete window offset
x1 = X(n) - floor(round(Uin(n))/2);
x2 = X(n) + ceil(round(Uin(n))/2);
y1 = Y(n) - floor(round(Vin(n))/2);
y2 = Y(n) + ceil(round(Vin(n))/2);
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);
for r=1:size(im1,3)
%find the image windows
zone1 = im1( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]),r);
zone2 = im2( max([1 ymin2]):min([L(1) ymax2]),max([1 xmin2]):min([L(2) xmax2]),r);
if size(zone1,1)~=Ny || size(zone1,2)~=Nx
w1 = zeros(Ny,Nx);
w1( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = zone1;
zone1 = w1;
end
if size(zone2,1)~=Ny || size(zone2,2)~=Nx
w2 = zeros(Ny,Nx);
w2( 1+max([0 1-ymin2]):Ny-max([0 ymax2-L(1)]),1+max([0 1-xmin2]):Nx-max([0 xmax2-L(2)]) ) = zone2;
zone2 = w2;
end
if Zeromean==1
zone1=zone1-mean(mean(zone1));
zone2=zone2-mean(mean(zone2));
end
%apply the image spatial filter
region1 = (zone1).*sfilt1;
region2 = (zone2).*sfilt2;
%FFTs and Cross-Correlation
f1 = fftn(region1,[Sy Sx]);
f2 = fftn(region2,[Sy Sx]);
P21 = f2.*conj(f1);
%Phase Correlation
W = ones(Sy,Sx);
Wden = sqrt(P21.*conj(P21));
W(Wden~=0) = Wden(Wden~=0);
if frac ~= 1
R = P21./(W.^frac); %Apply fractional weighting to the normalization
else
R = P21./W;
end
% If DRPC, the calculate the spectral function
% dynamically based on the autocorrelation
if dyn_rpc
CPS = ifftn(Wden,'symmetric');
[~,~,~,Drpc]=subpixel(CPS(fftindy,fftindx),Sx,Sy,cnorm,Peaklocator,0,D);
spectral = fftshift(energyfilt(Sx,Sy,Drpc/sqrt(2),0));
end
%Robust Phase Correlation with spectral energy filter
G = ifftn(R.*spectral,'symmetric');
G = G(fftindy,fftindx);
Gens(:,:,r) = abs(G);
end
G=mean(Gens,3);
%subpixel estimation
[U(n,:),V(n,:),Ctemp,Dtemp]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D);
if Peakswitch
C(n,:)=Ctemp;
Dia(n,:)=Dtemp;
end
if saveplane
Corrplanes(:,:,n) = G;
end
end
else
for n=1:length(X)
%apply the second order discrete window offset
x1 = X(n) - floor(round(Uin(n))/2);
x2 = X(n) + ceil(round(Uin(n))/2);
y1 = Y(n) - floor(round(Vin(n))/2);
y2 = Y(n) + ceil(round(Vin(n))/2);
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);
%find the image windows
zone1 = im1( max([1 ymin1]):min([L(1) ymax1]),max([1 xmin1]):min([L(2) xmax1]));
zone2 = im2( max([1 ymin2]):min([L(1) ymax2]),max([1 xmin2]):min([L(2) xmax2]));
if size(zone1,1)~=Ny || size(zone1,2)~=Nx
w1 = zeros(Ny,Nx);
w1( 1+max([0 1-ymin1]):Ny-max([0 ymax1-L(1)]),1+max([0 1-xmin1]):Nx-max([0 xmax1-L(2)]) ) = zone1;
zone1 = w1;
end
if size(zone2,1)~=Ny || size(zone2,2)~=Nx
w2 = zeros(Ny,Nx);
w2( 1+max([0 1-ymin2]):Ny-max([0 ymax2-L(1)]),1+max([0 1-xmin2]):Nx-max([0 xmax2-L(2)]) ) = zone2;
zone2 = w2;
end
if Zeromean==1
zone1=zone1-mean(mean(zone1));
zone2=zone2-mean(mean(zone2));
end
%apply the image spatial filter
region1 = (zone1).*sfilt1;
region2 = (zone2).*sfilt2;
%FFTs and Cross-Correlation
f1 = fftn(region1,[Sy Sx]);
f2 = fftn(region2,[Sy Sx]);
P21 = f2.*conj(f1);
%Phase Correlation
W = ones(Sy,Sx);
Wden = sqrt(P21.*conj(P21));
W(Wden~=0) = Wden(Wden~=0);
if frac ~= 1
R = P21./(W.^frac);%apply factional weighting to the normalization
else
R = P21./W;
end
% If DRPC, the calculate the spectral function
% dynamically based on the autocorrelation
if dyn_rpc
CPS = ifftn(Wden,'symmetric');
[~,~,~,Drpc]=subpixel(CPS(fftindy,fftindx),Sx,Sy,cnorm,Peaklocator,0,D);
spectral = fftshift(energyfilt(Sx,Sy,Drpc/sqrt(2),0));
end
%Robust Phase Correlation with spectral energy filter
G = ifftn(R.*spectral,'symmetric');
G = G(fftindy,fftindx);
G = abs(G);
%subpixel estimation
[U(n,:),V(n,:),Ctemp,Dtemp]=subpixel(G,Sx,Sy,cnorm,Peaklocator,Peakswitch,D);
if Peakswitch
C(n,:)=Ctemp;
Dia(n,:)=Dtemp;
end
if saveplane
Corrplanes(:,:,n) = G;
end
% Evaluate uncertainty options for RPC
if uncertainty.ppruncertainty(e)==1
%SNR calculation the other output arguments of Cal_SNR
%are Maximum peak value,PRMSR,PCE,ENTROPY
metric='PPR';
PPRval = Cal_SNR(G,metric);
% Save the SNR metrics
SNRmetric.PPR(n)=PPRval;
% Evaluate PPR Uncertainty
% John J Charonko Model
[Ux,Uy,~,~,~,~]=calibration_based_uncertainty('PPR_Charonkomodel',PPRval,upper(tcorr));
uncertainty2D.Upprx(n)=Ux;
uncertainty2D.Uppry(n)=Uy;
% Xue Model
[~,~,UxLB,UxUB,UyLB,UyUB]=calibration_based_uncertainty('PPR_Xuemodel',PPRval,upper(tcorr));
uncertainty2D.UpprxLB(n)=UxLB;
uncertainty2D.UppryLB(n)=UyLB;
uncertainty2D.UpprxUB(n)=UxUB;
uncertainty2D.UppryUB(n)=UyUB;
end
if uncertainty.miuncertainty(e)==1
%Autocorrelations
P11 = f1.*conj(f1);
P22 = f2.*conj(f2);
Auto1 = ifftn(P11,'symmetric');
Auto2 = ifftn(P22,'symmetric');
Auto1 = Auto1(fftindy,fftindx);
Auto2 = Auto2(fftindy,fftindx);
Auto1=abs(Auto1);
Auto2=abs(Auto2);
nAuto1 = Auto1-min(Auto1(:)); % Autocorrelation plane of image 1
nAuto2 = Auto2-min(Auto2(:)); % Autocorrelation plane of image 2
% 3 pt Gaussian fit to Autocorrelation Diameter
[~,~,~,~,Dauto1x3,Dauto1y3,~]=subpixel(nAuto1,Sx,Sy,cnorm,1,0,D);
[~,~,~,~,Dauto2x3,Dauto2y3,~]=subpixel(nAuto2,Sx,Sy,cnorm,1,0,D);
Diap1=sqrt(Dauto1x3*Dauto1y3/2);
Diap2=sqrt(Dauto2x3*Dauto2y3/2);
%Average Autocorrelation Diameter
Autod=mean([Diap1 Diap2]);
uncertainty2D.Autod(n)=Autod;
%MI Calculation
INTS1 = max(region1(:));
INTS2 = max(region2(:));
%Calculate the magnitude part of the correlation plane
W = ones(Sy,Sx);
Wden = sqrt(P21.*conj(P21));
W(Wden~=0) = Wden(Wden~=0);
% This part should be checked original function was a
% bit confusing but from the paper this should work
% W = ones(Sy,Sx);
% Wden = sqrt(P21.*conj(P21));
% Wden1 = ifftn(Wden,'symmetric');
% W(P21~=0) = Wden(P21~=0);
magG = ifftn(W,'symmetric');
magG = magG(fftindy,fftindx);
[MIval,~,~,~,~,~] = MI_Cal_RPC(magG,nAuto1,nAuto2,INTS1,INTS2,Dauto1x3,Dauto1y3,Dauto2x3,Dauto2y3,Sx,Sy,fftindx,fftindy);
SNRmetric.MI(n)=MIval;
% Estimate MI uncertainty
[~,~,UxLB,UxUB,UyLB,UyUB]=calibration_based_uncertainty('MI_Xuemodel',MIval,upper(tcorr));
uncertainty2D.UmixLB(n)=UxLB;
uncertainty2D.UmiyLB(n)=UyLB;
uncertainty2D.UmixUB(n)=UxUB;
uncertainty2D.UmiyUB(n)=UyUB;
end
% The uncertainty estimation using Moment of Correlation
% method for RPC needs to be figured out
% if uncertainty.mcuncertainty==1
% [MCx,MCy]=Moment_of_correlation(f1,f2,Sx,Sy,cnorm,D,fftindx,fftindy,G,DXtemp,DYtemp,region1,region2,Udiff,Vdiff);
% end
% % perform image matching uncertainty calculation (moved here by lalit)
% if uncertainty.imuncertainty(e)==1
% % Perform image matching uncertainty estimation
% [deltax,deltay,Npartu,Npartv] = original_particle_disparity(region1, region2, Sx/2, Sy/2, unique(res));
% % assign values to the uncertainty structure
% uncertainty2D.Uimx(n) = deltax;
% uncertainty2D.Uimy(n) = deltay;
% Npart(:,:,1) = Npartu;
% Npart(:,:,2) = Npartv;
% uncertainty2D.Nump(n) = mean(Npart,3);
% clear Npart
% end
% %% calculate uncertainty using correlation statistics (lalit)
% if uncertainty.csuncertainty(e) == 1
% [sigma_cs_x, sigma_cs_y] = correlation_statistics(region1, region2, window, res, Sx/2, Sy/2);
% % [sigma_cs_x, sigma_cs_y] = correlation_statistics(im1d, im2d, Wsize(e,:),Wres(:, :, e), X(Eval>=0),Y(Eval>=0));
% % [sigma_cs_x, sigma_cs_y] = correlation_statistics(im1, im2, Wsize(e,:),Wres(:, :, e), X(Eval>=0),Y(Eval>=0));
% % keyboard;
% % copy cs uncertainties in masked regions
% uncertainty2D.Ucsx(n) = sigma_cs_x;
% uncertainty2D.Ucsy(n) = sigma_cs_y;
% end
end
end
otherwise
%throw an error, we shouldn't be here
error('invalid correlation type')
end
%add DWO to estimation
U = round(Uin)+U;
V = round(Vin)+V;
end