diff --git a/Moment_of_correlation.m b/Moment_of_correlation.m new file mode 100644 index 0000000..6a4c4b4 --- /dev/null +++ b/Moment_of_correlation.m @@ -0,0 +1,236 @@ +function[Ixx,Iyy,biasx,biasy,Neff,Autod]=Moment_of_correlation(P21,f1,f2,Sx,Sy,cnorm,D,fftindx,fftindy,G,DXtemp,DYtemp,region1,region2,MIest) +%% This function calculates the standard uncertainty from the cross-correlation plane using Moment of Correlation method. +% written by Sayantan Bhattacharya + +% Inputs +%P21= cross-correlation plane in fourier domain +%f1= fft of window1/image1 +%f2= fft of window2/image2 +%Sx= window size Y +%Sy= window size X +%D= Approx Diameter of particle image (e.g. [2.8 2.8]), this is used to +%initialize the autocorrelation diameter estimation using subpixel fit + +%cnorm= Matrix of ones of size Sx,Sy this is non unity if you want to +%correct for spatial winfow filtering + +% fftindx,fftindy= fftshift indices I guess a simple fftshift can be used +% instead, but just to be consistent with prana functions + +%G= Cross correlation plane of two windows + +%region1, region2= spatially windowed image + +%DXtemp,DYtemp= Cross correlation peak diameter estimated from 3pt subpixel +%fit on cross-correlation plane G + +%MIest = If MI is already estimated then use that otherwise calculate MI + +%Output +%Ixx=PDF diameter in X direction +%Iyy=PDF diameter in X direction +%biasx=bias uncertainty in X direction +%biasy=bias uncertainty in Y direction +%Neff=Effective number of particles +%Autod= Average Autocorrelation diameter + + +if MIest==-1 + % If Mi has not been estimated + %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]); + + + %MI Calculation + INTS1 = max(region1(:)); + INTS2 = max(region2(:)); + [MI,~,~,~,~,~] = MI_Cal_SCC(G,nAuto1,nAuto2,INTS1,INTS2,Dauto1x3,Dauto1y3,Dauto2x3,Dauto2y3,Sx,Sy,fftindx,fftindy); +else + % Use estimated MI + MI=MIest; + Autod=0; +end +%Cross-correlation subpixel fit +% [U,V,~,~,DXtemp,DYtemp,~]=subpixel(G,Sx,Sy,cnorm,1,0,D); + + +if isnan(DXtemp) + DXtemp=sqrt(2)*2.8; +end +if isnan(DYtemp) + DYtemp=sqrt(2)*2.8; +end + +% Estimate the diameter of the cross-correlation peak (with initiation of +% 3pt fit estimates of correlation diameter) +[~,~,~,~,dxc1,dyc1,alpha2]=subpixel(G,Sx,Sy,cnorm,3,0,(1/sqrt(2))*[DXtemp DYtemp]); +% Taking care of peak rotation +DCCx = sqrt( (cos(alpha2)^2*dxc1^2 + sin(alpha2)^2*dyc1^2) ); +DCCy = sqrt( (sin(alpha2)^2*dxc1^2 + cos(alpha2)^2*dyc1^2) ); +% If estimates are too big than the 3pt fit estimate then default to 3 pt +% fit estimate of the correlation diameter +if DCCx >sqrt(2)*DXtemp + DCCx=DXtemp; +end +if DCCy >sqrt(2)*DYtemp + DCCy=DYtemp; +end + +%Use Cross-correlation diameter estimate to define convolving Gaussian Diameter +Dconv=(1/sqrt(2))*[DCCx DCCy]; + +% Finding the Phase correlation +W = ones(Sy,Sx); +Wden = sqrt(P21.*conj(P21)); +W(Wden~=0) = Wden(Wden~=0); +R = P21./W; % This is effectively the fourier transform of the PDF + + +% constructing the gaussian which will be convolved with the pdf +[Xt,Yt]=meshgrid(1:Sx,1:Sy); +gfilt=exp(-(4/Dconv(1)^2).*(Xt-Sx/2-1).^2-(4/Dconv(2)^2).*(Yt-Sy/2-1).^2); +Gfilt=abs(fftn(gfilt,[Sy Sx])); + +%Convolving the PDF and the Gaussian +G1 = ifftn(R.*Gfilt,'symmetric'); +G1 = G1(fftindy,fftindx); +G1 = abs(G1); +G1=G1-min(G1(:)); + + +%subpixel estimation using least squres guassfit for the convolved plane G1 +[gpx,gpy,~,~,dx1,dy1,alpha1]=subpixel(G1,Sx,Sy,cnorm,3,0,Dconv); + +%find the PDF major and minor axis +Px=real(((dx1(1)^2 - 2*(Dconv(1))^2)).^0.5); +Py=real(((dy1(1)^2 - 2*(Dconv(2))^2)).^0.5); + +%If the pdf diameter comes out to be zero or imaginary try 3point fit for +%the convolved plane G1 +if Px==0 || Py==0 + % if zero try 3point fit + [gpx,gpy,~,~,dx1,dy1,~]=subpixel(G1,Sx,Sy,cnorm,1,0,Dconv); + + Px=real(((dx1(1)^2 - 2*(Dconv(1))^2)).^0.5); + Py=real(((dy1(1)^2 - 2*(Dconv(2))^2)).^0.5); + alpha1=0; + %project major and minor axis to X and Y axes + Ixx = sqrt( 1/16 * (cos(alpha1)^2*Px^2 + sin(alpha1)^2*Py^2) ); + Iyy = sqrt( 1/16 * (sin(alpha1)^2*Px^2 + cos(alpha1)^2*Py^2) ); + + +else + % If not zero then continue + %project to axes + Ixx = sqrt( 1/16 * (cos(alpha1)^2*Px^2 + sin(alpha1)^2*Py^2) ); + Iyy = sqrt( 1/16 * (sin(alpha1)^2*Px^2 + cos(alpha1)^2*Py^2) ); + +end + +% The mean particle diameter is estimated from the mean cross-correlation +% peak diameter divided by sqrt(2) +part_dia=(1/sqrt(2))*mean([DCCx DCCy]); + +% Effective number of pixels is MI times a circular blob of pixels with +% mean particle diameter +Neff=MI*((pi/4)*(part_dia)^2); + +%Bias: This is the peak location of the convolved gaussian plane, typically +%this is done for multipass converged correlation plane in which case the +%peak should be at zero and any deviation is the bias. However this will +%not work for first pass. +biasx=gpx; +biasy=gpy; + +% The gradient correction and scaling is done outside PIVwindowed as the +% gradient estimation requires the full velcoity field + +% %Gradient correction using velocity gradient field Udiff and Vdiff +% Ixxt= real(sqrt(Ixx.^2 - (Autod.^2/16).*(Udiff).^2 )); +% Iyyt= real(sqrt(Iyy.^2 - (Autod.^2/16).*(Vdiff).^2 )); +% +% +% %MC uncertainty after scaling and bias correction +% MCx=sqrt(biasx^2+(Ixxt^2)/Neff); +% MCy=sqrt(biasy^2+(Iyyt^2)/Neff); +end + +% JJC: this function is an exact duplicate of the external function. +% Is it necessary? +function [MI,Nim1,Nim2,INTS,Diapx,Diapy] = MI_Cal_SCC(G,nAuto1,nAuto2,INTS1,INTS2,Diap1x,Diap1y,Diap2x,Diap2y,Sx,Sy,fftindx,fftindy) +%[MI,INTS,DiaP] = MI_Cal_SCC(G,nAuto1,nAuto2,INTS1,INTS2,Sx,Sy,fftindx,fftindy) +%[MI,INTS,DiaP] = MI_Cal_SCC( region1,region2,Sx,Sy,fftindx,fftindy) +%This function is used to calculate the mutual information between two +%consecutive frames using SCC method +% region 1 and 2 are the particle image within the interrogation window. +% Sx and Sy are the correlation window size +% fftindx and fftindy are fftshift indicies +% We will first find the maximum intensity from the image, and average +% particle diameter from particle image autocorrelation. A standard +% gaussian particle is built based on these two parameter. Then the +% contritbution of one particle is caluclate by take the auotcorrelatio +% peak of stardard gaussian particle. The primary peak on correlation +% plane is a summation of the contribution of all correlated particles. +% MI is the ratio between contribution of all correlated particles and +% one particle. + + + +% Diapx=(1/sqrt(2))*mean([Diap1x Diap2x]); +% Diapy=(1/sqrt(2))*mean([Diap1y Diap2y]); +Diapx=(1/sqrt(2))*(Diap1x*Diap2x)^0.5; +Diapy=(1/sqrt(2))*(Diap1y*Diap2y)^0.5; + +%make analytical gaussian particle and calculate +%autocorrelaiton of the standard particle +xco = 1:Sx;xco = repmat(xco,[Sx,1]); % build X axis +yco = 1:Sy;yco = repmat(yco',[1,Sy]); % build Y axis +% INTS = (INTS1+INTS2)/2; % intensity of standard Gaussian particle is the mean (maximum) intensity of two frames +INTS = sqrt(INTS1*INTS2); % intensity of standard Gaussian particle is the mean (maximum) intensity of two frames + +% % DiaP = (DiaP1+DiaP2)/2; % diameter of standard Gaussian particle is the mean diameter of two frames +% % fg = INTS*exp(-8*((xco-round(Sx/2)).^2+(yco-round(Sy/2)).^2)/DiaP^2); % build the particle intensity distribution based on 2D Gaussian distribution + +fg = INTS*exp(-8*((xco-round(Sx/2)).^2/(Diapx^2)+(yco-round(Sy/2)).^2/(Diapy^2))); % build the particle intensity distribution based on 2D Gaussian distribution + +fg = fftn(fg,[Sy,Sx]); %FFT +Pg = fg.*conj(fg); %FFT based correlation +Sp = ifftn(Pg,'symmetric'); % convert to time domain +Sp = Sp(fftindy,fftindx); % Sp is the autocorrelation of the standadrd Gaussian particle + + %use cross correaltion plane to get MI + Gnorm = G-min(G(:)); %minimum correlation subtraction to elminate background noise effect + [~,Gind] = max(Gnorm(:)); % find the primary peak location + shift_locyg = 1+mod(Gind-1,Sy); % find the Y coordinate of the peak + shift_locxg = ceil(Gind/Sy); % find the X coordinate of the peak + GXshift = Sy/2+1-shift_locxg; % find the distance between the peak location to the center of correlation plane in X direction + GYshift = Sx/2+1-shift_locyg; % find the distance between the peak location to the center of correlation plane in Y direction + SGnorm = circshift(Gnorm,[GYshift,GXshift]); % shift the peak to the center of the plane + %keyboard; +MI = max(SGnorm(:))/max(Sp(:)); % MI is the ratio between the contribution of all correlated particle and the contribution of one particle +Nim1=max(nAuto1(:))/max(Sp(:)); +Nim2=max(nAuto2(:))/max(Sp(:)); + + + +end + diff --git a/image_matching_prana.m b/image_matching_prana.m new file mode 100644 index 0000000..cdc1940 --- /dev/null +++ b/image_matching_prana.m @@ -0,0 +1,258 @@ +function [deltax,deltay,NPeak,dispx,dispy,cw]=image_matching_prana(region3,region4,dx,dy) +%THIS FUNCTION CALCULATES THE DISPARITY BETWEEN TWO CLOSELY MATCHING +%IMAGES(LIKE AFTER DWO OR ITERATIVE WINDOW DEFORMATION WHEN PARTICLES +%IN IMAGE PAIR IS CLOSE TO EACH OTHER). +%INPUT:- +%REGION3=WINDOW FROM IMAGE1, WITH ZERO MEAN AND 50% GAUSSIAN SPATIAL FILTER +%OPERATION APPLIED IN PRANA. +%REGION4=WINDOW FROM IMAGE2, WITH ZERO MEAN AND 50% GAUSSIAN SPATIAL FILTER +%OPERATION APPLIED IN PRANA. +%OUTPUT:- +%DELTAX=UNCERTAINTY IN X USING DISPARITY DISTRIBUTION DISPX +%DELTAY=UNCERTAINTY IN Y USING DISPARITY DISTRIBUTION DISPY +%THIS IS BASED ON THE PRINCIPLE DISCUSSED IN "PIV UNCERTAINTY QUANTIFICATION +%BY IMAGE MATCHING", MST. 2013.BY ANDREA SCIACCHITANO,BERNHARD WIENEKE AND +%FULVIO SCARANO +%THIS CODE WRITTEN FOR PRANA IMPLEMENTATION IS WRITTEN BY SAYANTAN BHATTACHARYA. + +if nargin<3 + dx = 0; + dy = 0; +end + +%% Threshold image and find peaks +sr=1; %SEARCH RADIUS FOR MATCHING PARTICLES +improduct=((region3).*(region4)); %PRODUCT OF INTENSITIES +rp3=improduct(improduct~=min(improduct(:)));%GETTING RID OF MINIMUM +thresh=0.5*rms(rp3(:));%SELECTING THRESHOLD 0.5 OF RMS OF INTENSITY PRODUCT + +improduct(improduct=(Ny-3) || Indx<=3 || Indx>=(Nx-3) + dispx(ji)=-100; + dispy(ji)=-100; + %cw(ji)=-100; + %fprintf('Limits exceeded'); + lflag='true'; + %keyboard; + continue; + end + + %FORM THE WINDOW ON WHICH TO SEARCH + sradiusy=Indy-sr:1:Indy+sr; + sradiusx=Indx-sr:1:Indx+sr; + + %SELECT THE REGION FROM INPUT IMAGES WHICH CORRESPONDS TO THE WINDOW + rA=region3(sradiusy,sradiusx); + rB=region4(sradiusy,sradiusx); + + %Find Peaks in regions rA and rB in 4 point neighbourhood + peakA=imregionalmax(rA,4);% 4 point neighbourhood + peakB=imregionalmax(rB,4);% 4 point neighbourhood + [im1y,im1x] =find(peakA==1); + [im2y,im2x] =find(peakB==1); + %PEAKS IN INDIVIDUAL IMAGE CORRESPONDING TO SEARCH WINDOW using maximum +% [im1y,im1x] =find(rA==max(rA(:))); +% [im2y,im2x] =find(rB==max(rB(:))); + + + %iF MULTIPLE PEAKS FOUND SELECT THE PEAK CLOSEST TO WINDOW CENTER + if length(im1y)>1 && length(im1x)>1 + %distance of peaks with respect to center + d1=((im1y-(sr+1)).^2+(im1x-(sr+1)).^2).^0.5; + in1=find(d1==min(d1));% find minimum + in1=in1(1); + %select minimum for window1 + im1y=im1y(in1); + im1x=im1x(in1); + + end + if length(im2y)>1 && length(im2x)>1 + %distance of peaks with respect to center + d2=((im2y-(sr+1)).^2+(im2x-(sr+1)).^2).^0.5; + in2=find(d2==min(d2));% find minimum + in2=in2(1); + %select minimum for window2 + im2y=im2y(in2); + im2x=im2x(in2); + end + + + %GETTING PEAK COORDINATE W.R.T. IMAGE COORDINATE SYSTEM + im1y=Indy+im1y-(sr+1); + im1x=Indx+im1x-(sr+1); + im2y=Indy+im2y-(sr+1); + im2x=Indx+im2x-(sr+1); + + %DOING SUBPIXEL FIT FOR EACH PEAK EACH DIRECTION + % for ith peak in image1, subpixel fit in y direction + lCm11 = log(region3(im1y-1,im1x)); + lC001 = log(region3(im1y,im1x)); + lCp11 = log(region3(im1y+1,im1x)); + if (2*(lCm11+lCp11-2*lC001)) == 0 % if three point Gaussian fit fail, set the particle diameter to a deafault guess value + shift_erry1 = 0; + else + shift_erry1 = (lCm11-lCp11)/(2*(lCm11+lCp11-2*lC001)); + end + clear lCm11 lC001 lCp11; + % for ith peak in image1, subpixel fit in x direction + lCm11 = log(region3(im1y,im1x-1)); + lC001 = log(region3(im1y,im1x)); + lCp11 = log(region3(im1y,im1x+1)); + if (2*(lCm11+lCp11-2*lC001)) == 0 % if three point Gaussian fit fail, set the particle diameter to a deafault guess value + shift_errx1 = 0; + else + shift_errx1 = (lCm11-lCp11)/(2*(lCm11+lCp11-2*lC001)); + end + clear lCm11 lC001 lCp11; + % for ith peak in image2, subpixel fit in y direction + lCm11 = log(region4(im2y-1,im2x)); + lC001 = log(region4(im2y,im2x)); + lCp11 = log(region4(im2y+1,im2x)); + if (2*(lCm11+lCp11-2*lC001)) == 0 % if three point Gaussian fit fail, set the particle diameter to a deafault guess value + shift_erry2 = 0; + else + shift_erry2 = (lCm11-lCp11)/(2*(lCm11+lCp11-2*lC001)); + end + clear lCm11 lC001 lCp11; + % for ith peak in image2, subpixel fit in x direction + lCm11 = log(region4(im2y,im2x-1)); + lC001 = log(region4(im2y,im2x)); + lCp11 = log(region4(im2y,im2x+1)); + if (2*(lCm11+lCp11-2*lC001)) == 0 % if three point Gaussian fit fail, set the particle diameter to a deafault guess value + shift_errx2 = 0; + else + shift_errx2 = (lCm11-lCp11)/(2*(lCm11+lCp11-2*lC001)); + end + clear lCm11 lC001 lCp11; + + + %PARTICLE CENTER LOCATIONS UPTO SUBPIXEL APPROXIMATION + Indy1=im1y+shift_erry1; + Indx1=im1x+shift_errx1; + Indy2=im2y+shift_erry2; + Indx2=im2x+shift_errx2; + + % If subpixel fit is nan or greater than 1 then subpixel fit failed, + % that point neglected + if abs(shift_errx1)>1 || abs(shift_errx2)>1 || abs(shift_erry1)>1 || abs(shift_erry2)>1 + dispx(ji)=-200; + dispy(ji)=-200; + elseif isnan(shift_errx1) || isnan(shift_errx2) || isnan(shift_erry1) || isnan(shift_erry2) + dispx(ji)=-300; + dispy(ji)=-300; + else + %DISPARITY BETWEEN PARTICLE PAIR FOR EACH PAIR OF MATCHED PARTICLE + dispx(ji)=Indx2-Indx1; + dispy(ji)=Indy2-Indy1; + end + +end + +%% DISCARDING POINTS WHICH WERE CLOSE TO THE BOUNDARY +if strcmp(lflag,'true') + %keyboard; + ind=find(dispx(:)==-100); + dispx(ind)=[]; + dispy(ind)=[]; + cw(ind)=[]; + %NPeak + NPeak=NPeak-length(ind); + %length(ind) + %NPeak +end +%% DISCARDING POINTS WHERE GAUSSIAN FIT FAILED +ind2=find(dispx(:)==-200); +dispx(ind2)=[]; +dispy(ind2)=[]; +cw(ind2)=[]; +%NPeak +NPeak=NPeak-length(ind2); + +ind3=find(dispy(:)==-300); +dispx(ind3)=[]; +dispy(ind3)=[]; +cw(ind3)=[]; +%NPeak +NPeak=NPeak-length(ind3); + +%% If after this elimination no Peak is detected return nan values for uncertainties +if NPeak==0 + deltax=nan; + deltay=nan; + return; +end + +%% Remove any subpixel shift left in the images +dispx = dispx - dx; +dispy = dispy - dy; + +%% CALCULATING WEIGHTED MEAN AND VARIANCE OF DISPARITY DISTRIBUTION FOR THIS +%WINDOW PAIRS AND SAVING THE UNCERTAINTY AS OUTPUT. +mewx=(1/sum(cw))*sum(cw.*dispx);% Weiighted mean +sigx=sqrt((1/sum(cw)).*sum(cw.*((dispx-mewx).^2)));% Weighted std deviation +deltax=sqrt(mewx^2 + (sigx/sqrt(NPeak))^2);% Uncertainty in x + +mewy=(1/sum(cw))*sum(cw.*dispy);% Weiighted mean +sigy=sqrt((1/sum(cw)).*sum(cw.*((dispy-mewy).^2)));% Weighted std deviation +deltay=sqrt(mewy^2 + (sigy/sqrt(NPeak))^2);% Uncertainty in y + +end \ No newline at end of file