diff --git a/stereo_uncertainty_codes/Calibration_job.mat b/stereo_uncertainty_codes/Calibration_job.mat
new file mode 100644
index 0000000..05be20d
Binary files /dev/null and b/stereo_uncertainty_codes/Calibration_job.mat differ
diff --git a/stereo_uncertainty_codes/IMcam13_ensemble_IM_uncertainty.mat b/stereo_uncertainty_codes/IMcam13_ensemble_IM_uncertainty.mat
new file mode 100644
index 0000000..a8ebbbc
Binary files /dev/null and b/stereo_uncertainty_codes/IMcam13_ensemble_IM_uncertainty.mat differ
diff --git a/stereo_uncertainty_codes/IMcam24_ensemble_IM_uncertainty.mat b/stereo_uncertainty_codes/IMcam24_ensemble_IM_uncertainty.mat
new file mode 100644
index 0000000..30ea0bb
Binary files /dev/null and b/stereo_uncertainty_codes/IMcam24_ensemble_IM_uncertainty.mat differ
diff --git a/stereo_uncertainty_codes/Plots/Thumbs.db b/stereo_uncertainty_codes/Plots/Thumbs.db
new file mode 100644
index 0000000..afbf144
Binary files /dev/null and b/stereo_uncertainty_codes/Plots/Thumbs.db differ
diff --git a/stereo_uncertainty_codes/Plots/Uhist.png b/stereo_uncertainty_codes/Plots/Uhist.png
new file mode 100644
index 0000000..b732fa0
Binary files /dev/null and b/stereo_uncertainty_codes/Plots/Uhist.png differ
diff --git a/stereo_uncertainty_codes/Plots/Vhist.png b/stereo_uncertainty_codes/Plots/Vhist.png
new file mode 100644
index 0000000..3bc61d2
Binary files /dev/null and b/stereo_uncertainty_codes/Plots/Vhist.png differ
diff --git a/stereo_uncertainty_codes/Plots/Whist.png b/stereo_uncertainty_codes/Plots/Whist.png
new file mode 100644
index 0000000..1cdf7de
Binary files /dev/null and b/stereo_uncertainty_codes/Plots/Whist.png differ
diff --git a/stereo_uncertainty_codes/Plots/desktop.ini b/stereo_uncertainty_codes/Plots/desktop.ini
new file mode 100644
index 0000000..64a4af2
--- /dev/null
+++ b/stereo_uncertainty_codes/Plots/desktop.ini
@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=This folder is shared online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=12
+
\ No newline at end of file
diff --git a/stereo_uncertainty_codes/Plots/planarhist.png b/stereo_uncertainty_codes/Plots/planarhist.png
new file mode 100644
index 0000000..f6c4431
Binary files /dev/null and b/stereo_uncertainty_codes/Plots/planarhist.png differ
diff --git a/stereo_uncertainty_codes/Read_me.docx b/stereo_uncertainty_codes/Read_me.docx
new file mode 100644
index 0000000..231dfc4
Binary files /dev/null and b/stereo_uncertainty_codes/Read_me.docx differ
diff --git a/stereo_uncertainty_codes/Thumbs.db b/stereo_uncertainty_codes/Thumbs.db
new file mode 100644
index 0000000..7309c3d
Binary files /dev/null and b/stereo_uncertainty_codes/Thumbs.db differ
diff --git a/stereo_uncertainty_codes/Uncertainty_propagation_through_mapping_function_gradient.m b/stereo_uncertainty_codes/Uncertainty_propagation_through_mapping_function_gradient.m
new file mode 100644
index 0000000..fd6cf31
--- /dev/null
+++ b/stereo_uncertainty_codes/Uncertainty_propagation_through_mapping_function_gradient.m
@@ -0,0 +1,149 @@
+function [Un_alpha1,Un_alpha2,Un_beta1,Un_beta2]=Uncertainty_propagation_through_mapping_function_gradient(calmat,xg,yg,zg,Uncalcoeff,unwx,unwy,unwz)
+% This function calculates the uncertainty in the mapping function gradients and subsequesntly the angle uncertainty
+%written by Sayantan Bhattacharya January 2016
+[r,c]=size(xg);
+dFdx1=zeros(r,c,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+dFdx2=zeros(r,c,4);
+dFdx3=zeros(r,c,4);
+Un_dFdx1=zeros(r,c,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+Un_dFdx2=zeros(r,c,4);
+Un_dFdx3=zeros(r,c,4);
+
+for gg=1:4
+ a=calmat(:,gg);
+ % Uc=uncoeff(:,gg);
+
+ % Mapping function gradient
+ dFdx1(:,:,gg) = a(2) + 2*a(5)*xg + a(6)*yg + a(8)*zg + 3*a(11)*xg.^2 + 2*a(12)*xg.*yg + ...
+ a(13)*yg.^2 + 2*a(15)*xg.*zg + a(16)*yg.*zg + a(18)*zg.^2;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*xg + 2*a(7)*yg + a(9)*zg + a(12)*xg.^2 + 2*a(13)*xg.*yg + ...
+ 3*a(14)*yg.^2 + a(16)*xg.*zg + 2*a(17)*yg.*zg + a(19)*zg.^2;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*xg + a(9)*yg + 2*a(10)*zg + a(15)*xg.^2 + a(16)*xg.*yg + ...
+ a(17)*yg.^2 + 2*a(18)*xg.*zg + 2*a(19)*yg.*zg;
+ % Second order derivative
+ cfxx=2*a(5) +6*a(11)*xg + 2*a(12)*yg + 2*a(15)*zg ;
+ cfxy=a(6) + 2*a(12)*xg + 2*a(13)*yg + a(16)*zg ;
+ cfxz=a(8)+ 2*a(15)*xg + a(16)*yg + 2*a(18)*zg;
+
+ cfyx=a(6) + 2*a(12)*xg + 2*a(13)*yg + a(16)*zg ;
+ cfyy=2*a(7) +2*a(13)*xg + 6*a(14)*yg + 2*a(17)*zg;
+ cfyz= a(9)+ a(16)*xg + 2*a(17)*yg + 2*a(19)*zg;
+
+ cfzx=a(8)+ 2*a(15)*xg + a(16)*yg + 2*a(18)*zg;
+ cfzy=a(9)+ a(16)*xg + 2*a(17)*yg + 2*a(19)*zg;
+ cfzz=2*a(10)+ 2*a(18)*xg + 2*a(19)*yg;
+
+ %Uncertainty due to Ux, Uy and Uz
+ T1dFdx1=sqrt((cfxx.^2).*(unwx.^2)+(cfxy.^2).*(unwy.^2)+(cfxz.^2).*(unwz.^2));
+ T1dFdx2=sqrt((cfyx.^2).*(unwx.^2)+(cfyy.^2).*(unwy.^2)+(cfyz.^2).*(unwz.^2));
+ T1dFdx3=sqrt((cfzx.^2).*(unwx.^2)+(cfzy.^2).*(unwy.^2)+(cfzz.^2).*(unwz.^2));
+
+
+ % Uncertainty due to cal coefficients
+
+ a=Uncalcoeff(:,gg); % Uncertainty in mapping function coefficients
+
+ T2dFdx1=sqrt(a(2).^2 + (2*a(5)).^2.*xg.^2 + a(6).^2*yg.^2 + a(8).^2*zg.^2 +...
+ (3*a(11)).^2*xg.^4 + (2*a(12)).^2*xg.^2.*yg.^2 + a(13).^2*yg.^4 + ...
+ (2*a(15)).^2*xg.^2.*zg.^2 + a(16).^2*yg.^2.*zg.^2 + a(18).^2*zg.^4);
+
+ T2dFdx2=sqrt(a(3).^2 + a(6).^2*xg.^2 + (2*a(7)).^2*yg.^2 + a(9).^2*zg.^2 +...
+ a(12).^2*xg.^4 + (2*a(13)).^2*xg.^2.*yg.^2 + (3*a(14)).^2*yg.^4 + ...
+ a(16).^2*xg.^2.*zg.^2 + (2*a(17)).^2*yg.^2.*zg.^2 + a(19).^2*zg.^4);
+
+ T2dFdx3=sqrt(a(4).^2 + a(8).^2*xg.^2 + a(9).^2*yg.^2 + (2*a(10)).^2*zg.^2 +...
+ a(15).^2*xg.^4 + a(16).^2*xg.^2.*yg.^2 + a(17).^2*yg.^4 + ...
+ (2*a(18)).^2*xg.^2.*zg.^2 + (2*a(19)).^2*yg.^2.*zg.^2);
+
+ %Total Uncertainty in mapping function gradient
+ Un_dFdx1(:,:,gg)=sqrt(T1dFdx1.^2+T2dFdx1.^2);
+ Un_dFdx2(:,:,gg)=sqrt(T1dFdx2.^2+T2dFdx2.^2);
+ Un_dFdx3(:,:,gg)=sqrt(T1dFdx3.^2+T2dFdx3.^2);
+
+ % Uncertainty in mapping function gradient due to Ux,Uy,Uz
+% Un_dFdx1(:,:,gg)=sqrt(T1dFdx1.^2);
+% Un_dFdx2(:,:,gg)=sqrt(T1dFdx2.^2);
+% Un_dFdx3(:,:,gg)=sqrt(T1dFdx3.^2);
+
+ % Uncertainty in mapping function gradient due to Uai's
+% Un_dFdx1(:,:,gg)=sqrt(T2dFdx1.^2);
+% Un_dFdx2(:,:,gg)=sqrt(T2dFdx2.^2);
+% Un_dFdx3(:,:,gg)=sqrt(T2dFdx3.^2);
+
+end
+% keyboard;
+
+%ANGLE CALCULATION
+N1=((dFdx3(:,:,2).*dFdx2(:,:,1))-(dFdx2(:,:,2).*dFdx3(:,:,1)));
+D1=(dFdx2(:,:,2).*dFdx1(:,:,1)-dFdx1(:,:,2).*dFdx2(:,:,1));
+alpha1 = N1./D1;
+%A1=atan2(N1,D1);
+
+N2=((dFdx3(:,:,4).*dFdx2(:,:,3))-(dFdx2(:,:,4).*dFdx3(:,:,3)));
+D2=(dFdx2(:,:,4).*dFdx1(:,:,3)-dFdx1(:,:,4).*dFdx2(:,:,3));
+alpha2 = N2./D2;
+%A2=atan2(N2,D2);
+
+N3=((dFdx3(:,:,2).*dFdx1(:,:,1))-(dFdx1(:,:,2).*dFdx3(:,:,1)));
+D3=(dFdx1(:,:,2).*dFdx2(:,:,1)-dFdx2(:,:,2).*dFdx1(:,:,1));
+beta1 = N3./D3;
+%B1=atan2(N3,D3);
+
+N4=((dFdx3(:,:,4).*dFdx1(:,:,3))-(dFdx1(:,:,4).*dFdx3(:,:,3)));
+D4=(dFdx1(:,:,4).*dFdx2(:,:,3)-dFdx2(:,:,4).*dFdx1(:,:,3));
+beta2 = N4./D4;
+%B2=atan2(N4,D4);
+
+
+%UNCERTAINTY IN ANGLES
+%uncertainty for alpha1
+mul1=(1./(1+alpha1.^2)).^2;
+term1=(((N1./(D1.^2)).*dFdx2(:,:,2)).^2).*(Un_dFdx1(:,:,1)).^2;
+term2=((dFdx3(:,:,2)./D1 + (N1./D1.^2).*dFdx1(:,:,2)).^2).*(Un_dFdx2(:,:,1)).^2;
+term3=((dFdx2(:,:,2)./D1).^2).*(Un_dFdx3(:,:,1)).^2;
+
+term4=(((N1./(D1.^2)).*dFdx2(:,:,1)).^2).*(Un_dFdx1(:,:,2)).^2;
+term5=((dFdx3(:,:,1)./D1 + (N1./D1.^2).*dFdx1(:,:,1)).^2).*(Un_dFdx2(:,:,2)).^2;
+term6=((dFdx2(:,:,1)./D1).^2).*(Un_dFdx3(:,:,2)).^2;
+
+Un_alpha1=(mul1.*(term1+term2+term3+term4+term5+term6)).^0.5;
+
+%uncertainty for alpha2
+mul2=(1./(1+alpha2.^2)).^2;
+term11=(((N2./(D2.^2)).*dFdx2(:,:,4)).^2).*(Un_dFdx1(:,:,3)).^2;
+term12=((dFdx3(:,:,4)./D2 + (N2./D2.^2).*dFdx1(:,:,4)).^2).*(Un_dFdx2(:,:,3)).^2;
+term13=((dFdx2(:,:,4)./D2).^2).*(Un_dFdx3(:,:,3)).^2;
+
+term14=(((N2./(D2.^2)).*dFdx2(:,:,3)).^2).*(Un_dFdx1(:,:,4)).^2;
+term15=((dFdx3(:,:,3)./D2 + (N2./D2.^2).*dFdx1(:,:,3)).^2).*(Un_dFdx2(:,:,4)).^2;
+term16=((dFdx2(:,:,3)./D2).^2).*(Un_dFdx3(:,:,4)).^2;
+
+Un_alpha2=(mul2.*(term11+term12+term13+term14+term15+term16)).^0.5;
+
+%uncertainty for beta1
+mul3=(1./(1+beta1.^2)).^2;
+term21=(((N3./(D3.^2)).*dFdx1(:,:,2)).^2).*(Un_dFdx2(:,:,1)).^2;
+term22=((dFdx3(:,:,2)./D3 + (N3./D3.^2).*dFdx2(:,:,2)).^2).*(Un_dFdx1(:,:,1)).^2;
+term23=((dFdx1(:,:,2)./D3).^2).*(Un_dFdx3(:,:,1)).^2;
+
+term24=(((N3./(D3.^2)).*dFdx1(:,:,1)).^2).*(Un_dFdx2(:,:,2)).^2;
+term25=((dFdx3(:,:,1)./D3 + (N3./D3.^2).*dFdx2(:,:,1)).^2).*(Un_dFdx1(:,:,2)).^2;
+term26=((dFdx1(:,:,1)./D3).^2).*(Un_dFdx3(:,:,2)).^2;
+
+Un_beta1=(mul3.*(term21+term22+term23+term24+term25+term26)).^0.5;
+
+%uncertainty for beta2
+mul4=(1./(1+beta2.^2)).^2;
+term131=(((N4./(D4.^2)).*dFdx1(:,:,4)).^2).*(Un_dFdx2(:,:,3)).^2;
+term132=((dFdx3(:,:,4)./D4 + (N4./D4.^2).*dFdx2(:,:,4)).^2).*(Un_dFdx1(:,:,3)).^2;
+term133=((dFdx1(:,:,4)./D4).^2).*(Un_dFdx3(:,:,3)).^2;
+
+term134=(((N4./(D4.^2)).*dFdx1(:,:,3)).^2).*(Un_dFdx2(:,:,4)).^2;
+term135=((dFdx3(:,:,3)./D4 + (N4./D4.^2).*dFdx2(:,:,3)).^2).*(Un_dFdx1(:,:,4)).^2;
+term136=((dFdx1(:,:,3)./D4).^2).*(Un_dFdx3(:,:,4)).^2;
+
+Un_beta2=(mul4.*(term131+term132+term133+term134+term135+term136)).^0.5;
+
+end
\ No newline at end of file
diff --git a/stereo_uncertainty_codes/VR13_scc_pass2_00025.mat b/stereo_uncertainty_codes/VR13_scc_pass2_00025.mat
new file mode 100644
index 0000000..40379e7
Binary files /dev/null and b/stereo_uncertainty_codes/VR13_scc_pass2_00025.mat differ
diff --git a/stereo_uncertainty_codes/Vortex_ring_stereo_uncertainty_using_CS.m b/stereo_uncertainty_codes/Vortex_ring_stereo_uncertainty_using_CS.m
new file mode 100644
index 0000000..ae90f98
--- /dev/null
+++ b/stereo_uncertainty_codes/Vortex_ring_stereo_uncertainty_using_CS.m
@@ -0,0 +1,464 @@
+%Uncertainty Propagation for Stereo Geometric Reconstruction(Willert
+%Method)for vortex ring, where planar uncertainties are evaluated using CS
+%method in Davis
+clc;
+clear;
+close all;
+
+
+lw=2;fs=18;
+
+printfig=0;
+
+main_dir=pwd;
+meth='CS';
+
+%Output Directory
+dirout=fullfile(main_dir,'Plots');
+if ~exist(dirout,'dir')
+ mkdir(dirout);
+end
+
+
+%Load Stereo calibration job
+stereojob= load(fullfile(main_dir,'Prana_stereo_job_after_selfcal.mat'));
+stdjob=stereojob.stdjob;
+%Calibration job
+caljob=stdjob.caljobfile.datasave.caljob;
+calmat=[caljob.aXcam1 caljob.aYcam1 caljob.aXcam2 caljob.aYcam2];
+calmat;
+%2D prana job
+planarjob=stdjob.pranajob.Data;
+
+%Load disparity field between two camera using ensemble correlation
+dispfield=load(fullfile(main_dir,'VR13_scc_pass2_00025.mat'));
+
+
+% Load Davis individual camera fields and stereo reconstructed fields
+Davissol= load(fullfile(main_dir,'Davis_processed_result.mat'));
+Davis13=Davissol.Davis13;
+
+%Load true solution
+truesol=load(fullfile(main_dir,'Vortex_Ring_True_Solution_interp_on_solution_grid.mat'));
+Ut=truesol.Ut;
+Vt=truesol.Vt;
+Wt=truesol.Wt;
+%Magnification in mm/pix
+scaling=stdjob.scaling.wil;
+mx=scaling.xscale;
+my=scaling.yscale;
+% Here the results are for camera1camera3 pair
+cam1='1';cam2='3';
+if strcmp(cam1,'1') && strcmp(cam2,'3')
+ mz=scaling.yscale;
+elseif strcmp(cam1,'2') && strcmp(cam2,'4')
+ mz=scaling.xscale;
+end
+
+%seconds per frame
+spf=0.001;
+
+%No. of frames
+t=50;
+
+%Starting frame
+fstart=24;
+
+
+%% Find uncertainty in angles
+
+%directory which has uncertainty in disparity field
+disp_uncertainty_filedir=main_dir;
+
+[Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,tanalpha1,tanalpha2,tanbeta1,tanbeta2]=stereo_angle_uncertainty(caljob,planarjob,dispfield,disp_uncertainty_filedir);
+
+%%
+[Sx,Sy]=size(tanalpha1);
+
+%Initialization
+Us=zeros(Sx,Sy,t);Vs=zeros(Sx,Sy,t);Ws=zeros(Sx,Sy,t);
+errx1=zeros(Sx,Sy,t);erry1=zeros(Sx,Sy,t);errz1=zeros(Sx,Sy,t);
+imx1=zeros(Sx,Sy,t);imy1=zeros(Sx,Sy,t);
+imx2=zeros(Sx,Sy,t);imy2=zeros(Sx,Sy,t);
+unim_u=zeros(Sx,Sy,t);unim_v=zeros(Sx,Sy,t);unim_w=zeros(Sx,Sy,t);
+
+for i=1:t
+
+ %% Load reconstructed solution
+ Us(:,:,i)=(spf/mx)*1e3*Davis13.Uw(1:end-1,1:end-1,i);
+ Vs(:,:,i)=(spf/my)*1e3*Davis13.Vw(1:end-1,1:end-1,i);
+ Ws(:,:,i)=(spf/mz)*1e3*Davis13.Ww(1:end-1,1:end-1,i);
+
+ xgrid=1000*Davis13.Xw(1:end-1,1:end-1);
+ ygrid=1000*Davis13.Yw(1:end-1,1:end-1);
+ zgrid=1000*Davis13.Zw(1:end-1,1:end-1);
+
+% [r,c]=size(xgrid);
+
+
+ %% Get error fields
+ errx1(:,:,i)=(Ut(:,:,i)-Us(:,:,i));
+ erry1(:,:,i)=(Vt(:,:,i)-Vs(:,:,i));
+ errz1(:,:,i)=(Wt(:,:,i)-Ws(:,:,i));
+
+
+
+ %% Load cam1cam2 2d fields
+ %get 2d velocities
+ U1=Davis13.Us1(1:end-1,1:end-1,i);
+ V1=Davis13.Vs1(1:end-1,1:end-1,i);
+ X1=Davis13.X(1:end-1,1:end-1);
+ Y1=Davis13.Y(1:end-1,1:end-1);
+
+ U2=Davis13.Us2(1:end-1,1:end-1,i);
+ V2=Davis13.Vs2(1:end-1,1:end-1,i);
+ X2=Davis13.X(1:end-1,1:end-1);
+ Y2=Davis13.Y(1:end-1,1:end-1);
+
+ %Get 2d uncertainties
+ %Camera1
+ imx1(:,:,i)=Davis13.UCSx1(1:end-1,1:end-1,i);%;
+ imy1(:,:,i)=Davis13.UCSy1(1:end-1,1:end-1,i);%
+
+ %Camera2
+ imx2(:,:,i)=Davis13.UCSx2(1:end-1,1:end-1,i);%
+ imy2(:,:,i)=Davis13.UCSy2(1:end-1,1:end-1,i);%
+
+
+ %% Stereo Uncertainty Propagation
+
+ Unu1=imx1(:,:,i);
+ Unu2=imx2(:,:,i);
+ Unv1=imy1(:,:,i);
+ Unv2=imy2(:,:,i);
+
+ [Un_u1,Un_v1,Un_w1,JU,JV,JW]=stereo_uncertainty_propagation(Unu1,Unv1,Unu2,Unv2,U1,V1,U2,V2,Ws(:,:,i),Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,tanalpha1,tanalpha2,tanbeta1,tanbeta2,mx,my);
+
+
+ unim_u(:,:,i)=Un_u1;
+ unim_v(:,:,i)=Un_v1;
+ unim_w(:,:,i)=Un_w1;
+
+end
+
+%Good measurement cut off
+ulx=0.5;uly=0.5;ulz=1.5;
+%% Plot RMS Uncertainty and Error fields (Spatial RMS Profile)
+
+for i=1:Sx
+ for j=1:Sy
+
+
+ rms2derrx=squeeze(abs(errx1(i,j,:)));
+ rms2dunx=squeeze(abs(unim_u(i,j,:)));
+ ix1=find(rms2derrx>ulx);
+ rms2derrx(ix1)=[];
+ rms2dunx(ix1)=[];
+% Nx(i,j)=length(ix1);
+
+ rms2derry=squeeze(abs(erry1(i,j,:)));
+ rms2duny=squeeze(abs(unim_v(i,j,:)));
+ iy1=find(rms2derry>uly);
+ rms2derry(iy1)=[];
+ rms2duny(iy1)=[];
+% Ny(i,j)=length(iy1);
+
+ rms2derrz=squeeze(abs(errz1(i,j,:)));
+ rms2dunz=squeeze(abs(unim_w(i,j,:)));
+ iz1=find(rms2derrz>ulz);
+ rms2derrz(iz1)=[];
+ rms2dunz(iz1)=[];
+% Nz(i,j)=length(iz1);
+
+
+ Ue(i,j)=rms(rms2derrx);
+ Uun(i,j)=rms(rms2dunx);
+ Ve(i,j)=rms(rms2derry);
+ Vun(i,j)=rms(rms2duny);
+ We(i,j)=rms(rms2derrz);
+ Wun(i,j)=rms(rms2dunz);
+
+ end
+end
+
+figure;
+subplot(2,3,1);imagesc(Ue);caxis([0 0.2]);
+subplot(2,3,2);imagesc(Ve);caxis([0 0.2]);
+subplot(2,3,3);imagesc(We);caxis([0 0.45]);
+subplot(2,3,4);imagesc(Uun);caxis([0 0.2]);
+subplot(2,3,5);imagesc(Vun);caxis([0 0.2]);
+subplot(2,3,6);imagesc(Wun);caxis([0 0.45]);
+%% Plot RMS Uncertainty and Error fields (Temporal RMS Profile)
+
+for k=1:t
+
+ rex=squeeze(abs(errx1(:,:,k)));
+ rIMx=squeeze(abs(unim_u(:,:,k)));
+ rex=rex(:);
+ rIMx=rIMx(:);
+ ix1=find(rex>ulx);
+ rex(ix1)=[];
+ rIMx(ix1)=[];
+% Nx(i,j)=length(ix1);
+
+ rey=squeeze(abs(erry1(:,:,k)));
+ rIMy=squeeze(abs(unim_v(:,:,k)));
+ rey=rey(:);
+ rIMy=rIMy(:);
+ iy1=find(rey>uly);
+ rey(iy1)=[];
+ rIMy(iy1)=[];
+% Ny(i,j)=length(iy1);
+
+ rez=squeeze(abs(errz1(:,:,k)));
+ rIMz=squeeze(abs(unim_w(:,:,k)));
+ rez=rez(:);
+ rIMz=rIMz(:);
+ iz1=find(rez>ulz);
+ rez(iz1)=[];
+ rIMz(iz1)=[];
+% Nz(i,j)=length(iz1);
+
+
+ siguerr(k)=rms(rex);
+ siguun(k)=rms(rIMx);
+
+ sigverr(k)=rms(rey);
+ sigvun(k)=rms(rIMy);
+
+ sigwerr(k)=rms(rez);
+ sigwun(k)=rms(rIMz);
+
+
+end
+
+figure;hold on;
+plot(1:t,siguerr,'r-o',1:t,siguun,'r-*');
+plot(1:t,sigverr,'b-o',1:t,sigvun,'b-*');
+plot(1:t,sigwerr,'g-o',1:t,sigwun,'g-*');
+hold off;
+title('RMS Time Series (Cam24)');xlabel('Frame No.');
+ylabel('RMS error, uncertainty(pix/frame)');
+% hleg=legend({'\sigma^{e}_{u}','\sigma^{CS}_{u}','\sigma^{e}_{v}','\sigma^{CS}_{v}','\sigma^{e}_{w}','\sigma^{CS}_{w}'});
+% set(hleg,'location','eastoutside');
+set(gca,'FontSize',20);
+axis([0 50 0 0.5]);axis square;
+%% RMS 2D Uncertainty spatial maps
+nc=20;
+IMX1=rms(imx1,3);
+IMY1=rms(imy1,3);
+IMX2=rms(imx2,3);
+IMY2=rms(imy2,3);
+figure;
+subplot(2,2,1);contour(xgrid,ygrid,IMX1,nc);caxis([0 0.45]);title('\sigma^{CS}_{U1}');axis square;
+subplot(2,2,2);contour(xgrid,ygrid,IMY1,nc);caxis([0 0.45]);title('\sigma^{CS}_{V1}');axis square;
+subplot(2,2,3);contour(xgrid,ygrid,IMX2,nc);caxis([0 0.45]);title('\sigma^{CS}_{U2}');axis square;
+subplot(2,2,4);contour(xgrid,ygrid,IMY2,nc);caxis([0 0.45]);title('\sigma^{CS}_{V2}');axis square;
+
+%%
+% save([meth,'_Cam',cam1,cam2,'_RMSprofiles.mat'],'xgrid','ygrid','IMX1',...
+% 'IMX2','IMY1','IMY2','siguerr','siguun','sigverr','sigvun','sigwerr','sigwun',...
+% 'Ue','Uun','Ve','Vun','We','Wun');
+
+%%
+%[mean(errx(:)) mean(erry(:)) mean(errz(:))]
+
+% [mean(imx1(:)) mean(imy1(:)) mean(Nim1(:)) mean(imx2(:)) mean(imy2(:)) mean(Nim2(:))]
+% [rms(imx1(:)) rms(imy1(:)) mean(Nim1(:)) rms(imx2(:)) rms(imy2(:)) mean(Nim2(:))]
+
+%%
+errx=(errx1);
+erry=(erry1);
+errz=(errz1);
+
+% errx1(abs(errx1)>0.3)=[];
+% erry1(abs(erry1)>0.3)=[];
+% errz1(abs(errz1)>0.3)=[];
+%
+% Nbin1=60;
+% figure;hist(errx1(:),Nbin1);
+% figure;hist(erry1(:),Nbin1);
+% figure;hist(errz1(:),Nbin1);
+
+Nc=100;
+
+ul=0.5;
+
+vx1=linspace(0,ul,Nc);
+vy1=linspace(0,ul,Nc);
+
+vx2=linspace(0,ul,Nc);
+vy2=linspace(0,ul,Nc);
+
+Nu1=histc(imx1(:),vx1);
+Nv1=histc(imy1(:),vy1);
+Nu2=histc(imx2(:),vx2);
+Nv2=histc(imy2(:),vy2);
+
+figure;set(gcf,'DefaultLineLineWidth',1);set(gca,'FontSize',fs);
+plot(vx1,Nu1,'r-o',vy1,Nv1,'b-o',vx2,Nu2,'r-+',vy2,Nv2,'b-+');
+title('Planar Correlation Uncertainties');
+legend('U1','V1','U2','V2');
+if printfig==1
+ print(gcf,'-dpng',fullfile(dirout,'planarhist.png'),'-r300');
+end
+%% eliminate invalid measurements
+errx=errx(:);
+erry=erry(:);
+errz=errz(:);
+unim_u=unim_u(:);
+unim_v=unim_v(:);
+unim_w=unim_w(:);
+
+
+% %Good measurement cut off
+% ulx=0.5;uly=0.5;ulz=1.5;
+
+inx=find(abs(errx(:))>ulx);
+iny=find(abs(erry(:))>uly);
+inz=find(abs(errz(:))>ulz);
+
+errx(inx)=[];
+erry(iny)=[];
+errz(inz)=[];
+unim_u(inx)=[];
+unim_v(iny)=[];
+unim_w(inz)=[];
+
+
+
+% [mean(errx(:)) mean(erry(:)) mean(errz(:))]
+%% Calculate coverage
+
+cnt1=0;cnt2=0;cnt3=0;
+l1=length(errx);
+l2=length(erry);
+l3=length(errz);
+for m1=1:l1
+ if abs(errx(m1))<=unim_u(m1)
+ cnt1=cnt1+1;
+ end
+
+end
+
+for m2=1:l2
+ if abs(erry(m2))<=unim_v(m2)
+ cnt2=cnt2+1;
+ end
+
+end
+
+for m3=1:l3
+
+ if abs(errz(m3))<=unim_w(m3)
+ cnt3=cnt3+1;
+ end
+
+end
+covx1=100*cnt1/l1;
+covy1=100*cnt2/l2;
+covz1=100*cnt3/l3;
+fprintf('\n Coverage with angle uncertainty \n Ucov=%2.2f Vcov=%2.2f Wcov=%2.2f \n',[covx1 covy1 covz1]);
+%% Plot Histograms
+% if ispc
+% addpath Z:\Planar_Uncertainty_work\codes\export_fig\;
+% addpath Z:\Planar_Uncertainty_work\codes\color_code\;
+% else
+% addpath /home/shannon/a/bhattac3/Planar_Uncertainty_work/codes/export_fig/;
+% addpath /home/shannon/a/bhattac3/Planar_Uncertainty_work/codes/color_code/;
+% end
+
+set(0,'DefaultAxesFontName', 'Times New Roman');
+
+%Plot histograms
+Nb=30;
+
+vecx=linspace(0,ulx,Nb);
+vecy=linspace(0,uly,Nb);
+vecz=linspace(0,ulz,Nb);
+vecx1=linspace(-ulx,ulx,2*Nb-1);
+vecy1=linspace(-uly,uly,2*Nb-1);
+vecz1=linspace(-ulz,ulz,2*Nb-1);
+
+Nex=histc(errx,vecx1)/length(errx);
+Nimx=histc(unim_u,vecx)/length(errx);
+
+Ney=histc(erry,vecy1)/length(erry);
+Nimy=histc(unim_v,vecy)/length(erry);
+
+Nez=histc(errz,vecz1)/length(errz);
+Nimz=histc(unim_w,vecz)/length(errz);
+
+
+% X uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecx1,Nex,'-k');
+plot(vecx,Nimx,'r-');
+% HL1=fill(vecx1,Nex,'-k');
+% HL2=fill(vecx,Nimx,'r-');
+% set(HL1,'facealpha',.5);
+% set(HL2,'facealpha',.5);
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigu=rms(errx(abs(errx)<=ulx)).*ones(length(ll),1);
+mixx1=rms(unim_u(unim_u<=ulx)).*ones(length(ll),1);
+plot(sigu,ll,'k--',mixx1,ll,'r--');
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{x}',['\sigma_{x}^{',meth,'}']);
+axis([-ulx ulx 0 inf]);
+% title(['U Uncertainty Histogram coverage= ',num2str(covx1,'%02.02f')]);
+title('U Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+ print(gcf,'-dpng',fullfile(dirout,'Uhist.png'),'-r300');
+% export_fig(gcf,fullfile([dirout,'Uhist.png']),'-painters','-r360');
+end
+
+% Y uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecy1,Ney,'-k');
+plot(vecy,Nimy,'b-');
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigv=rms(erry(abs(erry)<=uly)).*ones(length(ll),1);
+miyy1=rms(unim_v(unim_v<=uly)).*ones(length(ll),1);
+plot(sigv,ll,'k--',miyy1,ll,'b--');
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{y}',['\sigma_{y}^{',meth,'}']);
+axis([-uly uly 0 inf]);
+% title(['V uncertainty histogram coverage= ',num2str(covy1,'%02.02f')]);
+title('V Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+ print(gcf,'-dpng',fullfile(dirout,'Vhist.png'),'-r300');
+% export_fig(gcf,fullfile([dirout,'Vhist.png']),'-painters','-r360');
+end
+
+% Z uncertainty Histogram
+figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
+plot(vecz1,Nez,'-k');
+plot(vecz,Nimz,'g-');
+yl=ylim(gca);
+ll=linspace(0,max(yl(2)),10);
+sigw=rms(errz(abs(errz)<=ulz)).*ones(length(ll),1);
+mizz1=rms(unim_w(unim_w<=ulz)).*ones(length(ll),1);
+plot(sigw,ll,'k--',mizz1,ll,'g--');
+
+grid on;box on; set(gcf,'color','white');hold off;
+
+legend('e_{z}',['\sigma_{z}^{',meth,'}']);
+axis([-ulz ulz 0 inf]);
+% title(['W uncertainty histogram coverage= ',num2str(covz1,'%02.02f')]);
+title('W Uncertainty Histogram');
+xlabel('RMS error and Uncertainty (pixels)');ylabel('No. of vectors');
+if printfig==1
+ print(gcf,'-dpng',fullfile(dirout,'Whist.png'),'-r300');
+% export_fig(gcf,fullfile([dirout,'Whist.png']),'-painters','-r360');
+end
+fprintf('\n RMS error & Uncertainty for x,y,z');
+fprintf('\n ex,ux,ey,uy,ez,uz \n %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f \n',[sigu(1) mixx1(1) sigv(1) miyy1(1) sigw(1) mizz1(1)]);
+
diff --git a/stereo_uncertainty_codes/desktop.ini b/stereo_uncertainty_codes/desktop.ini
new file mode 100644
index 0000000..64a4af2
--- /dev/null
+++ b/stereo_uncertainty_codes/desktop.ini
@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=This folder is shared online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=12
+
\ No newline at end of file
diff --git a/stereo_uncertainty_codes/imagedewarp_alone.m b/stereo_uncertainty_codes/imagedewarp_alone.m
new file mode 100644
index 0000000..4dfbade
--- /dev/null
+++ b/stereo_uncertainty_codes/imagedewarp_alone.m
@@ -0,0 +1,581 @@
+function [outputdirlist,dewarp_grid,scaling]=imagedewarp_alone(caldata,dewarpmethod,imagelist,imagesize,vectorlist)
+%This code dewarps the images or vector grid depending on the
+%reconstruction type and outputs the dewarped common grid coordinates and
+%the modified magnification or scaling
+%Inputs:-
+%caldata=structure containing all the calibration information
+%dewarpmethod='Soloff' or 'Willert
+%imagelist=list containing directories of individual camera images (for Willert or selfcal)
+%vectorlist=list containing directories of individual camera vector fields(Soloff)
+%Outputs:-
+%outputdirlist=dewarped image output directories
+%dewarp_grid=dewarped grid coordinates
+%scaling=magnification based on dewarped grid
+
+% This file is part of prana, an open-source GUI-driven program for
+% calculating velocity fields using PIV or PTV.
+%
+% Copyright (C) 2014 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 .
+
+%keyboard;
+orderz=caldata.modeltype;
+alldata.orderz=orderz;
+optionsls=caldata.optionsls;
+
+aXcam1=caldata.aXcam1;
+aYcam1=caldata.aYcam1;
+aXcam2=caldata.aXcam2;
+aYcam2=caldata.aYcam2;
+
+%keyboard;
+
+%JJC: We will be dewarping images?
+if strcmp(dewarpmethod,'Willert')
+
+ %keyboard;
+ dir1=imagelist.imdirec;
+ dir2=imagelist.imdirec2;
+ base1=imagelist.imbase;
+ base2=imagelist.imbase2;
+ ext=imagelist.imext;
+ zer=str2double(imagelist.imzeros);
+% fstep=str2double(imagelist.imfstep);
+ fstart=str2double(imagelist.imfstart);
+ fend=str2double(imagelist.imfend);
+ cstep=str2double(imagelist.imcstep);
+ dirsave=imagelist.outdirec;
+ dirout1=fullfile(dirsave,['Dewarped Images1',filesep]);
+ dirout2=fullfile(dirsave,['Dewarped Images2',filesep]);
+
+ if ~exist(dirout1,'dir')
+ mkdir(dirout1);
+ end
+ if ~exist(dirout2,'dir')
+ mkdir(dirout2);
+ end
+
+
+ outputdirlist.dewarpdir1=dirout1;
+ outputdirlist.dewarpdir2=dirout2;
+
+ istring1=sprintf(['%%s%%0%0.0fd.',ext],zer);
+ %keyboard;
+% IML=imread(fullfile(dir1,sprintf(istring1,base1,fstart)));
+% IMR=imread(fullfile(dir2,sprintf(istring1,base2,fstart)));
+% [Jmax1,Imax1] = size(IML); %Possible bug, isn't this [NY,NX], but used below as [NX,NY] to build X1points, etc.?
+% [Jmax2,Imax2] = size(IMR);
+ %bug fixed?
+ Imax1= imagesize(1); Jmax1= imagesize(2);
+ Imax2= imagesize(1); Jmax2= imagesize(2);
+
+ %The corners of each image?
+ %[SW,SE,NW,NE] ?
+
+ % this assigns the corner points of the image depending on the
+ % targetside i.e. if both the cameras are on the same side or opposite side
+ if caldata.targetsidecam1==caldata.targetsidecam2
+ X1points=[1 Jmax1 1 Jmax1];
+ Y1points=[1 1 Imax1 Imax1];
+ X2points=[1 Jmax2 1 Jmax2];
+ Y2points=[1 1 Imax2 Imax2];
+ elseif caldata.targetsidecam1==1 && caldata.targetsidecam2==2
+ X1points=[1 Jmax1 1 Jmax1];
+ Y1points=[1 1 Imax1 Imax1];
+ X2points=[Jmax2 1 Jmax2 1];
+ Y2points=[1 1 Imax2 Imax2];
+ elseif caldata.targetsidecam1==2 && caldata.targetsidecam2==1
+ X1points=[Jmax1 1 Jmax1 1];
+ Y1points=[1 1 Imax1 Imax1];
+ X2points=[1 Jmax2 1 Jmax2];
+ Y2points=[1 1 Imax2 Imax2];
+ end
+
+ % X1points=[0.5 Jmax1-0.5 0.5 Jmax1-0.5];
+ % Y1points=[0.5 0.5 Imax1-0.5 Imax1-0.5];
+ % X2points=[0.5 Jmax2-0.5 0.5 Jmax2-0.5];
+ % Y2points=[0.5 0.5 Imax2-0.5 Imax2-0.5];
+
+%JJC: We will be dewarping vector fields?
+elseif strcmp(dewarpmethod,'Soloff')
+ veclist1 = load(vectorlist{1});
+ x1=veclist1.X;
+ y1=veclist1.Y;
+ clear veclist1;
+
+ veclist2 = load(vectorlist{2});
+ x2=veclist2.X;
+ y2=veclist2.Y;
+ clear veclist2;
+
+ % %Possible bug, isn't this [NY,NX], but used below as [NX,NY] to build X1points, etc.?
+ % [Jmax1, Imax1]=size(x1);
+ % [Jmax2, Imax2]=size(x2);
+ %bug fixed?
+ [Imax1, Jmax1]=size(x1);
+ [Imax2, Jmax2]=size(x2);
+
+ outputdirlist='';
+
+ %Caution: These positions come from vector fields, and therefore may be
+ % vector coordinates that don't match the pixel coordinates!
+ % Do we need to translate them?
+ %To make this work (ignoring coordinate system mismatch for now), we
+ % must be assuming that the vector field locations are in units of
+ % pixels in Image space, not physical units and not physical (object?)
+ % space.
+ %
+ %The corners of each vector field?
+ %[SW,SE,NW,NE]? x is flipped E-W for backside cameras
+ % this assigns the corner points of the vector grid depending on the
+ % targetside i.e. if both the cameras are on the same side or opposite side
+ if caldata.targetsidecam1==caldata.targetsidecam2
+ X1points=[min(min(x1)) max(max(x1)) min(min(x1)) max(max(x1))] + 0.5;
+ Y1points=[min(min(y1)) min(min(y1)) max(max(y1)) max(max(y1))] + 0.5;
+ X2points=[min(min(x2)) max(max(x2)) min(min(x2)) max(max(x2))] + 0.5;
+ Y2points=[min(min(y2)) min(min(y2)) max(max(y2)) max(max(y2))] + 0.5;
+ % different arrangement if camera 2 is on other side
+ elseif caldata.targetsidecam1==1 && caldata.targetsidecam2==2
+
+ X1points=[min(min(x1)) max(max(x1)) min(min(x1)) max(max(x1))] + 0.5;
+ Y1points=[min(min(y1)) min(min(y1)) max(max(y1)) max(max(y1))] + 0.5;
+ X2points=[max(max(x2)) min(min(x2)) max(max(x2)) min(min(x2))] + 0.5;
+ Y2points=[min(min(y2)) min(min(y2)) max(max(y2)) max(max(y2))] + 0.5;
+
+ elseif caldata.targetsidecam1==2 && caldata.targetsidecam2==1
+
+ X1points=[max(max(x1)) min(min(x1)) max(max(x1)) min(min(x1))] + 0.5;
+ Y1points=[min(min(y1)) min(min(y1)) max(max(y1)) max(max(y1))] + 0.5;
+ X2points=[min(min(x2)) max(max(x2)) min(min(x2)) max(max(x2))] + 0.5;
+ Y2points=[min(min(y2)) min(min(y2)) max(max(y2)) max(max(y2))] + 0.5;
+
+ end
+
+end
+
+
+x0=[1 1]; % initial guess for solver
+% finds the over lap for the two cameras. This is only performed the
+% first time and then the information is used subsequently.
+%if c==1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Determination of Area Overlap
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% find overlapping area using the corner points of the first loaded
+% images, the order for X1points and similar is [bl br tl tr]
+
+fprintf('Computing overlapping area...\n');
+
+bottom1X=linspace(X1points(1),X1points(2),Jmax1)';
+top1X=linspace(X1points(3),X1points(4),Jmax1)';
+left1X=X1points(1)*ones(1,Imax1)';
+right1X=X1points(2)*ones(1,Imax1)';
+bottom1Y=Y1points(1)*ones(1,Jmax1)';
+top1Y=Y1points(3)*ones(1,Jmax1)';
+left1Y=linspace(Y1points(1),Y1points(3),Imax1)';
+right1Y=linspace(Y1points(2),Y1points(4),Imax1)';
+
+bottom2X=linspace(X2points(1),X2points(2),Jmax2)';
+top2X=linspace(X2points(3),X2points(4),Jmax2)';
+left2X=X2points(1)*ones(1,Imax2)';
+right2X=X2points(2)*ones(1,Imax2)';
+bottom2Y=Y2points(1)*ones(1,Jmax2)';
+top2Y=Y2points(3)*ones(1,Jmax2)';
+left2Y=linspace(Y2points(1),Y2points(3),Imax2)';
+right2Y=linspace(Y2points(2),Y2points(4),Imax2)';
+
+%keyboard;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% calculate x,y for the bottom,top,left,right vectors for X,Y from
+% camera 1
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+alldata.aX=aXcam1';
+alldata.aY=aYcam1';
+
+bottom1xy = zeros(2,Jmax1);
+top1xy = zeros(2,Jmax1);
+points=[bottom1X bottom1Y];
+for k=1:Jmax1
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 1
+ [bottom1xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+points=[top1X top1Y];
+for k=1:Jmax1
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 1
+ [top1xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+
+left1xy = zeros(2,Imax1);
+right1xy = zeros(2,Imax1);
+points=[left1X left1Y];
+for k=1:Imax1
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 1
+ [left1xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+points=[right1X right1Y];
+for k=1:Imax1
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 1
+ [right1xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% calculate x,y for the bottom,top,left,right vectors for X,Y from
+% camera 2
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+alldata.aX=aXcam2';
+alldata.aY=aYcam2';
+
+bottom2xy = zeros(2,Jmax2);
+top2xy = zeros(2,Jmax2);
+points=[bottom2X bottom2Y];
+for k=1:Jmax2
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 2
+ [bottom2xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+points=[top2X top2Y];
+for k=1:Jmax2
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 2
+ [top2xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+
+left2xy = zeros(2,Imax2);
+right2xy = zeros(2,Imax2);
+points=[left2X left2Y];
+for k=1:Imax2
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 2
+ [left2xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+points=[right2X right2Y];
+for k=1:Imax2
+ alldata.XYpoint=points(k,:)';
+ % solve for x,y for camera 2
+ [right2xy(:,k),~,~]=fsolve(@(x) poly_3xy_123z_2eqns(x,alldata),x0,optionsls);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% make object coordinate grid
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+xlow=max([left1xy(1,:) left2xy(1,:)]);
+xhigh=min([right1xy(1,:) right2xy(1,:)]);
+ylow=max([bottom1xy(2,:) bottom2xy(2,:)]);
+yhigh=min([top1xy(2,:) top2xy(2,:)]);
+
+%set the number of points in dewarped image to the size of Camera 1's
+%original undewarped image.
+ImaxD = Imax1; %number of points in y
+JmaxD = Jmax1; %number of points in x
+
+% [xgrid,ygrid] = meshgrid(linspace(-55,55,1024),linspace(-50,50,1024));
+[xgrid,ygrid]=meshgrid(linspace(xlow,xhigh,JmaxD),linspace(ylow,yhigh,ImaxD));
+%zgrid=zeros(size(xgrid));
+overplots=1;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Plot fig to check overlap
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+skp=32;
+figure;
+if overplots == 1
+ % use this to plot overlapping area in x,y
+
+ H(1) = plot(bottom1xy(1,:),bottom1xy(2,:),'.');hold on;
+ plot(top1xy(1,:),top1xy(2,:),'.');hold on;
+ plot(left1xy(1,:),left1xy(2,:),'.');hold on;
+ plot(right1xy(1,:),right1xy(2,:),'.');hold on;
+ H(2) = plot(bottom2xy(1,:),bottom2xy(2,:),'.');hold on;
+ plot(top2xy(1,:),top2xy(2,:),'.');hold on;
+ plot(left2xy(1,:),left2xy(2,:),'.');hold on;
+ plot(right2xy(1,:),right2xy(2,:),'.');hold on;
+ H(3) = plot([xlow xlow xhigh xhigh xlow],[ylow yhigh yhigh ylow ylow],'-k','LineWidth',1);hold on;
+ H2 = plot(xgrid(1:skp:end,1:skp:end),ygrid(1:skp:end,1:skp:end),'.r','MarkerSize',4);
+ xlabel('x (mm)');ylabel('y (mm)');
+ H(4) = H2(1);
+ title('Camera Overlap and new vector locations');
+ Lstr = {'Camera 1 border','Camera 2 border','Overlap Border','Vector location'};
+ legend(H,Lstr);
+ % set(L,'Position',[0.4 0.4 0.2314 0.1869])
+ % slashlocs = find(data.outputdirectory == '/');
+ % set(gcf,'name',data.outputdirectory(slashlocs(end)+1:end))
+ % axis([-200 100 -150 250])
+ drawnow
+end
+hold off;
+% print(gcf,'-dpng','GridOverlap.png','-r300');
+%angle calculation
+% [rows,cols]=size(xgrid);
+%
+% xgrid=X1;
+% ygrid=Y1;
+% zgrid=zeros(size(xgrid));
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % Compute gradients of calibration functions
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% aall=[caldata.aXcam1 caldata.aYcam1 caldata.aXcam2 caldata.aYcam2];
+% dFdx1=zeros(rows,cols,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+% dFdx2=zeros(rows,cols,4);
+% dFdx3=zeros(rows,cols,4);
+%
+%
+% if caldata.modeltype==1
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % Mapping the camera coord. to the World Coord. using 1sr order z
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% for gg=1:4
+% a=aall(:,gg);
+% dFdx1(:,:,gg) = a(2) + 2*a(5)*xgrid + a(6)*ygrid + a(8)*zgrid + 3*a(10)*xgrid.^2 + ...
+% 2*a(11)*xgrid.*ygrid + a(12)*ygrid.^2 + 2*a(14)*xgrid.*zgrid + a(15)*ygrid.*zgrid;
+%
+% dFdx2(:,:,gg) = a(3) + a(6)*xgrid + 2*a(7)*ygrid + a(9)*zgrid + a(11)*xgrid.^2 + ...
+% 2*a(12)*xgrid.*ygrid + 3*a(13)*ygrid.^2 + a(15)*xgrid.*zgrid + 2*a(16)*ygrid.*zgrid;
+%
+% dFdx3(:,:,gg) = a(4) + a(8)*xgrid + a(9)*ygrid + a(14)*xgrid.^2 + a(15)*xgrid.*ygrid + a(16)*ygrid.^2;
+% end
+% elseif caldata.modeltype==2
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% % Mapping the camera coord. to the World Coord. using 2nd order z
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% for gg=1:4
+% a=aall(:,gg);
+% dFdx1(:,:,gg) = a(2) + 2*a(5)*xgrid + a(6)*ygrid + a(8)*zgrid + 3*a(11)*xgrid.^2 + 2*a(12)*xgrid.*ygrid + ...
+% a(13)*ygrid.^2 + 2*a(15)*xgrid.*zgrid + a(16)*ygrid.*zgrid + a(18)*zgrid.^2;
+%
+% dFdx2(:,:,gg) = a(3) + a(6)*xgrid + 2*a(7)*ygrid + a(9)*zgrid + a(12)*xgrid.^2 + 2*a(13)*xgrid.*ygrid + ...
+% 3*a(14)*ygrid.^2 + a(16)*xgrid.*zgrid + 2*a(17)*ygrid.*zgrid + a(19)*zgrid.^2;
+%
+% dFdx3(:,:,gg) = a(4) + a(8)*xgrid + a(9)*ygrid + 2*a(10)*zgrid + a(15)*xgrid.^2 + a(16)*xgrid.*ygrid + ...
+% a(17)*ygrid.^2 + 2*a(18)*xgrid.*zgrid + 2*a(19)*ygrid.*zgrid;
+% end
+% end
+%
+% alpha1 = ((dFdx3(:,:,2).*dFdx2(:,:,1))-(dFdx2(:,:,2).*dFdx3(:,:,1)))./(dFdx2(:,:,2).*dFdx1(:,:,1)-dFdx1(:,:,2).*dFdx2(:,:,1));
+% beta1 = ((dFdx3(:,:,2).*dFdx1(:,:,1))-(dFdx1(:,:,2).*dFdx3(:,:,1)))./(dFdx1(:,:,2).*dFdx2(:,:,1)-dFdx2(:,:,2).*dFdx1(:,:,1));
+%
+% alpha2 = ((dFdx3(:,:,4).*dFdx2(:,:,3))-(dFdx2(:,:,4).*dFdx3(:,:,3)))./(dFdx2(:,:,4).*dFdx1(:,:,3)-dFdx1(:,:,4).*dFdx2(:,:,3));
+% beta2 = ((dFdx3(:,:,4).*dFdx1(:,:,3))-(dFdx1(:,:,4).*dFdx3(:,:,3)))./(dFdx1(:,:,4).*dFdx2(:,:,3)-dFdx2(:,:,4).*dFdx1(:,:,3));
+
+%Display camera angles for reference
+
+% figure(100); subplot(2,2,1);
+% imagesc(atand(alpha1)); colorbar; %caxis([25 30]);
+% title('Camera 1 Angle \alpha1','FontSize',16)
+% subplot(2,2,2);
+% imagesc(atand(alpha2)); colorbar; %caxis([-30 -25]);
+% title('Camera 2 Angle \alpha2','FontSize',16)
+% subplot(2,2,3);
+% imagesc(atand(beta1)); colorbar; %caxis([-2 2]);(end-(zed+10):end-(zed+5))
+% title('Camera 1 Angle \beta1','FontSize',16)
+% subplot(2,2,4);
+% imagesc(atand(beta2)); colorbar; %caxis([-2 2]);
+% title('Camera 1 Angle \beta1','FontSize',16)
+
+% [mean(atand(alpha1(:))) mean(atand(alpha2(:))) mean(atand(beta1(:))) mean(atand(beta2(:)))]
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Compute this grid in the IMAGE (object) plane to interpolate values
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[Xgrid1,Ygrid1]=poly_3xy_123z_fun(xgrid,ygrid,orderz,aXcam1,aYcam1);
+[Xgrid2,Ygrid2]=poly_3xy_123z_fun(xgrid,ygrid,orderz,aXcam2,aYcam2);
+
+dewarp_grid.Xgrid1=Xgrid1;
+dewarp_grid.Ygrid1=Ygrid1;
+dewarp_grid.Xgrid2=Xgrid2;
+dewarp_grid.Ygrid2=Ygrid2;
+dewarp_grid.xgrid=xgrid;
+dewarp_grid.ygrid=ygrid;
+
+% keyboard;
+%scaling is only an estimate returned as an output, and is not used in the
+% dewarping calculation
+%divide extent in object space by number of pixels.
+if strcmp(dewarpmethod,'Willert')
+ %JJC: this looks like a bug. Commented out and replaced with fix below
+ % % Does this have an off-by-one error? If we have four pixels at 1:4 with
+ % % unit scaling, max_x-min_x = 4-1 = 3, and Imax = 4, then we get
+ % % xscale = 3/4 = 0.75, not 1.0 like we'd expect.
+ % scaling.xscale =(max(xgrid(:))-min(xgrid(:)))/Imax1;
+ % scaling.yscale =(max(ygrid(:))-min(ygrid(:)))/Jmax1;
+
+ %Both numerator and denominator need to be (Max - Min), see Soloff
+ %below for analogous construction
+ scaling.xscale =(max(xgrid(:))-min(xgrid(:)))/(Jmax1-1);
+ scaling.yscale =(max(ygrid(:))-min(ygrid(:)))/(Imax1-1);
+elseif strcmp(dewarpmethod,'Soloff')
+ %JJC: this was already correct, see above
+ scaling.xscale =(max(xgrid(:))-min(xgrid(:)))/(max(x1(:))-min(x1(:)));
+ scaling.yscale =(max(ygrid(:))-min(ygrid(:)))/(max(y1(:))-min(y1(:)));
+end
+
+%keyboard;
+if strcmp(dewarpmethod,'Willert')
+ %[x1,y1] = meshgrid(0.5:1:Jmax1-0.5,0.5:1:Imax1-0.5);
+ dewarpflag=0; % added dewarpflag which is by default true but one can make it 0 to not run it if required
+ if dewarpflag==1
+ fprintf('Dewarping Images...\n');
+ delete(gcp);
+ parpool; % calling matlabpool for dewarping in parallel but in future versions it should be parpool
+ parfor k=fstart:fend+1
+ %reading recorded images
+
+ IMLi= im2double(imread(fullfile(dir1,sprintf(istring1,base1,k))));
+ IMRi= im2double(imread(fullfile(dir2,sprintf(istring1,base2,k+cstep-1))));
+
+ incl=class(IMLi);
+ %flipping images
+ IMLi=IMLi(end:-1:1,:);
+ IMRi=IMRi(end:-1:1,:);
+ %Interpolating on a common grid
+ %sincBlackmanInterp2 assumes images are on grid [1 NX] and [1 NY]
+% IMLo=((sincBlackmanInterp2(IMLi, Xgrid1, Ygrid1, 8,'blackman')));
+% IMRo=((sincBlackmanInterp2(IMRi, Xgrid2, Ygrid2, 8,'blackman')));
+
+ IMLo=whittaker_blackman(IMLi, Xgrid1, Ygrid1, 6,1);
+ IMRo=whittaker_blackman(IMRi, Xgrid2, Ygrid2, 6,1);
+
+
+
+ % IMLo=double((interp2(x1,y1,IMLi, Xgrid1, Ygrid1, 'spline',0)));
+ % IMRo=double((interp2(x1,y1,IMRi, Xgrid2, Ygrid2, 'spline',0)));
+
+
+ % keyboard;
+ %flipping them back for saving
+ IMLo=IMLo(end:-1:1,:);
+ IMRo=IMRo(end:-1:1,:);
+ IMLo=cast(IMLo,incl);
+ IMRo=cast(IMRo,incl);
+
+ % figure(4);imagesc(IMLo);colormap('gray');%axis equal tight;
+ % figure(5);imagesc(IMRo);colormap('gray');%axis equal tight;
+ % figure(6);imagesc(IMLi);colormap('gray');%axis equal tight;
+ % figure(7);imagesc(IMRi);colormap('gray');%axis equal tight;
+ %writing dewarped images to the output directory
+ %writing the images because of two reasons:-
+ %1)prana assumes all images to be in a sigle directory while processing
+ %2)In prana processing just the image matrices cannot be given as input, have to give image location
+ % clear IM;
+ % IM(:,:,1)=IMLo;
+ % IM(:,:,2)=IMRo;
+ % IM(:,:,3)=IMRo;
+ % figure(10);imshow(IM);
+
+ %keyboard;
+ %JJC: why the heck are we saving them with no TIFF compression?
+ imwrite((IMLo),fullfile(dirout1,sprintf(istring1,base1,k)),'TIFF','WriteMode','overwrite','Compression','none');
+ imwrite((IMRo),fullfile(dirout2,sprintf(istring1,base2,k+cstep-1)),'TIFF','WriteMode','overwrite','Compression','none');
+ %keyboard;
+
+ end
+ delete(gcp);
+ end
+end
+
+end
+
+function F=poly_3xy_123z_2eqns(x,alldata)
+% F=poly_3xy_123z_2eqns(x,alldata)
+% this function solves for the xy object coordinates with input
+% image coordiantes alldata.XYpoint. the resulting x vector contains
+% the (x y) object coords. This is for S-PIV so the z coord. is 0.
+
+% This function is called by reconstructvectorsfun.m
+
+% Writen by M. Brady
+% Edited and Commented by S. Raben
+
+aX=alldata.aX;
+aY=alldata.aY;
+orderz=alldata.orderz;
+XYpoint=alldata.XYpoint;
+
+if orderz==1 % cubic xy, linear z
+ polylist=[1 x(1) x(2) 0 x(1)^2 x(1)*x(2) x(2)^2 0 0 x(1)^3 x(1)^2*x(2) x(1)*x(2)^2 x(2)^3 0 0 0]';
+ Fpoly=[aX*polylist;aY*polylist]-XYpoint;
+
+elseif orderz==2 % cubic xy, quadratic z
+ polylist=[1 x(1) x(2) 0 x(1)^2 x(1)*x(2) x(2)^2 0 0 0 x(1)^3 x(1)^2*x(2) x(1)*x(2)^2 x(2)^3 0 0 0 0 0]';
+ Fpoly=[aX*polylist;aY*polylist]-XYpoint;
+
+else % camera pinhole
+
+end
+
+F=Fpoly;
+end
+
+function [Xgrid,Ygrid]=poly_3xy_123z_fun(xgrid,ygrid,orderz,aX,aY)
+% [Xgrid Ygrid]=poly_3xy_123z_fun(xgrid,ygrid,orderz,aX,aY)
+%
+
+% Writen by M. Brady
+% Edited and Commented by S. Raben
+
+x1=xgrid;
+x2=ygrid;
+[r,c]=size(xgrid);
+x3=zeros(r,c);
+
+if orderz==1 % cubic xy, linear z
+
+ Xgrid=aX(1) + aX(2).*x1 + aX(3).*x2 + aX(4).*x3 + aX(5).*x1.^2 +...
+ aX(6).*x1.*x2 + aX(7).*x2.^2 + aX(8).*x1.*x3 + aX(9).*x2.*x3 +...
+ aX(10).*x1.^3 + aX(11).*x1.^2.*x2 + aX(12).*x1.*x2.^2 +...
+ aX(13).*x2.^3 + aX(14).*x1.^2.*x3 + aX(15).*x1.*x2.*x3 +...
+ aX(16).*x2.^2.*x3;
+
+ Ygrid=aY(1) + aY(2).*x1 + aY(3).*x2 + aY(4).*x3 + aY(5).*x1.^2 +...
+ aY(6).*x1.*x2 + aY(7).*x2.^2 + aY(8).*x1.*x3 + aY(9).*x2.*x3 +...
+ aY(10).*x1.^3 + aY(11).*x1.^2.*x2 + aY(12).*x1.*x2.^2 +...
+ aY(13).*x2.^3 + aY(14).*x1.^2.*x3 + aY(15).*x1.*x2.*x3 +...
+ aY(16).*x2.^2.*x3;
+
+elseif orderz==2 % cubic xy, quadratic z
+
+ Xgrid=aX(1) + aX(2).*x1 + aX(3).*x2 + aX(4).*x3 + aX(5).*x1.^2 +...
+ aX(6).*x1.*x2 + aX(7).*x2.^2 + aX(8).*x1.*x3 + aX(9).*x2.*x3 +...
+ aX(10).*x3.^2 + aX(11).*x1.^3 + aX(12).*x1.^2.*x2 + aX(13).*x1.*x2.^2 +...
+ aX(14).*x2.^3 + aX(15).*x1.^2.*x3 + aX(16).*x1.*x2.*x3 +...
+ aX(17).*x2.^2.*x3 + aX(18).*x1.*x3.^2 + aX(19).*x2.*x3.^2;
+
+ Ygrid=aY(1) + aY(2).*x1 + aY(3).*x2 + aY(4).*x3 + aY(5).*x1.^2 +...
+ aY(6).*x1.*x2 + aY(7).*x2.^2 + aY(8).*x1.*x3 + aY(9).*x2.*x3 +...
+ aY(10).*x3.^2 + aY(11).*x1.^3 + aY(12).*x1.^2.*x2 + aY(13).*x1.*x2.^2 +...
+ aY(14).*x2.^3 + aY(15).*x1.^2.*x3 + aY(16).*x1.*x2.*x3 +...
+ aY(17).*x2.^2.*x3 + aY(18).*x1.*x3.^2 + aY(19).*x2.*x3.^2;
+
+else % pinhole
+
+
+end
+end
+
diff --git a/stereo_uncertainty_codes/stereo_angle_uncertainty.m b/stereo_uncertainty_codes/stereo_angle_uncertainty.m
new file mode 100644
index 0000000..d9b3a1a
--- /dev/null
+++ b/stereo_uncertainty_codes/stereo_angle_uncertainty.m
@@ -0,0 +1,310 @@
+function [Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,tanalpha1,tanalpha2,tanbeta1,tanbeta2]=stereo_angle_uncertainty(caljob,planarjob,dispfield,disp_uncertainty_filedir)
+%This function calculates the stereo angle uncertainties
+
+% Input Variables
+
+%caljob= This is the calibration job structure generated through Prana
+%Calibration and conatins the calibration mapping function coefficients and
+%the calibration parameters.
+%caljob.aXcam1 caljob.aYcam1 caljob.aXcam2 caljob.aYcam2 gives the X and Y
+%mapping function coefficents (polynomial model) for camera 1 and 2. This
+%can be before or after self-calibration and can be obtained using Prana.
+
+%caljob.modeltype=1 for polynomial model (cubic in x and y and linear in z)
+%caljob.modeltype=2 for polynomial model (cubic in x and y and quadratic in z)
+
+%world coordinates (from prana calibration)
+%caljob.allx1data= [x,y,z] world coordinates of calibration grid points for camera1
+%caljob.allx2data= [x,y,z] world coordinates of calibration grid points for camera2
+%caljob.allX1data= [X,Y] image coordinates of calibration grid points for camera1
+%caljob.allX2data= [X,Y] image coordinates of calibration grid points for camera2
+
+%planarjob= The 2D prana jobfile for correlating each camera image pair,
+%This argument contains the individual camera imagelist which is used to dewarp the images using the function
+%"imagedewarp_alone()". This input can be replaced by "dewarp_grid"
+%variable, which essentially contains the dewarped x and y grid on the
+%common region of overlap between two cameras.
+
+%dispfield is a structure which contains the disparity field between
+%two cameras. The fields in the structure are U,V,X and Y corresponding to
+%a disparity vector map obtained using ensemble correlation of two camera
+%images at the same time instant.
+
+
+%disp_uncertainty_filedir= The directory to store the calculated uncertainty in each disparity vector or
+%load the uncertainty from a previously calculated matfile in this directory.
+
+% Output Variables
+% tanalpha1,tanalpha2,tanbeta1,tanbeta2= tangent of the stereo angles alpha
+% (angle in x-z plane) and beta (angle in y-z plane) for camera 1 and 2
+
+% Un_alpha1,Un_alpha2,Un_beta1,Un_beta2= the uncertainty in stereo
+% angles alpha1, beta1, alpha2, beta2
+
+% written by Sayantan Bhattacharya January 2016
+
+%% Get Existing Calibration Matrix
+%This assigns the existing calibration matrix to calmat variable
+calmat=[caljob.aXcam1 caljob.aYcam1 caljob.aXcam2 caljob.aYcam2];
+
+%% Uncertainty estimated using ensmeble image matching
+% Here you load the estimated uncertainty for each disparity vector.
+% This is obtained in an ensemble sense i.e. first the disparity field for
+% the existing calibration is obtianed using ensemble correlation of camera
+% dewarped image pairs. The disparity field is then used to shift the individual
+% camera images onto each other(similar to image matching technique) and
+% obtain an error distribution for each particle position. The errors are
+% stacked up in a histogram for all particles in a window and over all
+% frames for which the disparity vector is evaluated. From the ensemble
+% error distribution of particle position mismatch the uncertainty in x and
+% y for each disparity vector is determined. This is in essence an ensemble
+% image matching technique.
+
+loadensmdspfld=load(fullfile(disp_uncertainty_filedir,'IMcam13_ensemble_IM_uncertainty.mat'));
+Unwx=loadensmdspfld.Unwx;% Uncertainty in X-direction for disparity vector
+Unwy=loadensmdspfld.Unwy;% Uncertainty in Y-direction for disparity vector
+dispfield.Unwx=Unwx;
+dispfield.Unwy=Unwy;
+
+%% Get Uncertainty in calibration coefficients and the world coordinates xyz
+%% using uncertainty propagation through Triangulation (Both bias and random uncertainty)
+%Get dewarped x and y coordinate grid on the overlapping camera domain
+imsize=[1024 1024];
+[~,dewarp_grid,~]=imagedewarp_alone(caljob,'Willert',planarjob,imsize);
+xgrid1=dewarp_grid.xgrid;
+ygrid1=dewarp_grid.ygrid;
+
+%Get z projected locations based on both mean and random uncertainties
+%evaluated with ensemble processing
+[xg,yg,Uncalcoeff,Ux,Uy,Uz,~,~]=uncertainty_in_xyz_calcoeff(caljob,dispfield,xgrid1,ygrid1);
+
+zg=zeros(size(xg));
+
+%% calculate the uncertainty in the angles using uncertainty in x,y,x and calibration coefficients and the mapping function gradient
+
+% calculate angles
+[tanalpha1,tanalpha2,tanbeta1,tanbeta2]=calculate_angle(calmat,xg,yg,zg,caljob.modeltype);
+
+% Find uncertainty in mapping fuction gradient and stereo angles
+[Un_alpha1,Un_alpha2,Un_beta1,Un_beta2]=Uncertainty_propagation_through_mapping_function_gradient(calmat,xg,yg,zg,Uncalcoeff,Ux,Uy,Uz);
+
+% figure;
+% subplot(2,2,1);imagesc(rad2deg(atan(alpha1)));colorbar;title('alpha1');
+% subplot(2,2,2);imagesc(rad2deg(atan(alpha2)));colorbar;title('alpha2');
+% subplot(2,2,3);imagesc(rad2deg(atan(beta1)));colorbar;title('beta1');
+% subplot(2,2,4);imagesc(rad2deg(atan(beta2)));colorbar;title('beta2');
+% print(gcf,'-dpng',[outputdir,'angles.png'],'-r300');
+
+% figure;
+% subplot(2,2,1);imagesc(rad2deg(Un_alpha1));colorbar;title('Un_{alpha1}');caxis([0 0.8]);
+% subplot(2,2,2);imagesc(rad2deg(Un_alpha2));colorbar;title('Un_{alpha2}');caxis([0 0.8]);
+% subplot(2,2,3);imagesc(rad2deg(Un_beta1));colorbar;title('Un_{beta1}');caxis([0 0.8]);
+% subplot(2,2,4);imagesc(rad2deg(Un_beta2));colorbar;title('Un_{beta2}');caxis([0 0.8]);
+% print(gcf,'-dpng',[outputdir,'angle_uncertainty.png'],'-r300');
+
+
+
+end
+function [tanalpha1,tanalpha2,tanbeta1,tanbeta2]=calculate_angle(calmat,xg,yg,zg,modeltype)
+% This function calculates the tangent of the camera angles based on the
+% mapping function gradient
+%written by Sayantan Bhattacharya January 2016
+% % % Viewing angles for camera 1 calculated at gridpoints xgrid,ygrid;
+[rows,cols]=size(xg);
+
+aall=calmat;%[aXcam1 aYcam1 aXcam2 aYcam2];
+dFdx1=zeros(rows,cols,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+dFdx2=zeros(rows,cols,4);
+dFdx3=zeros(rows,cols,4);
+
+if modeltype==1
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Mapping the camera coord. to the World Coord. using 1sr order z
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for gg=1:4
+ a=aall(:,gg);
+ dFdx1(:,:,gg) = a(2) + 2*a(5)*xg + a(6)*yg + a(8)*zg + 3*a(10)*xg.^2 + ...
+ 2*a(11)*xg.*yg + a(12)*yg.^2 + 2*a(14)*xg.*zg + a(15)*yg.*zg;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*xg + 2*a(7)*yg + a(9)*zg + a(11)*xg.^2 + ...
+ 2*a(12)*xg.*yg + 3*a(13)*yg.^2 + a(15)*xg.*zg + 2*a(16)*yg.*zg;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*xg + a(9)*yg + a(14)*xg.^2 + a(15)*xg.*yg + a(16)*yg.^2;
+ end
+
+elseif modeltype==2
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Mapping the camera coord. to the World Coord. using 2nd order z
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for gg=1:4
+ a=aall(:,gg);
+ dFdx1(:,:,gg) = a(2) + 2*a(5).*xg + a(6)*yg + a(8)*zg + 3*a(11)*xg.^2 + 2*a(12)*xg.*yg + ...
+ a(13)*yg.^2 + 2*a(15)*xg.*zg + a(16)*yg.*zg + a(18)*zg.^2;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*xg + 2*a(7)*yg + a(9)*zg + a(12)*xg.^2 + 2*a(13)*xg.*yg + ...
+ 3*a(14)*yg.^2 + a(16)*xg.*zg + 2*a(17)*yg.*zg + a(19)*zg.^2;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*xg + a(9)*yg + 2*a(10)*zg + a(15)*xg.^2 + a(16)*xg.*yg + ...
+ a(17)*yg.^2 + 2*a(18)*xg.*zg + 2*a(19)*yg.*zg;
+ end
+
+end
+
+
+tanalpha1=(((dFdx3(:,:,2).*dFdx2(:,:,1)) - (dFdx2(:,:,2).*dFdx3(:,:,1)))./((dFdx2(:,:,2).*dFdx1(:,:,1)) - (dFdx1(:,:,2).*dFdx2(:,:,1))));
+tanbeta1=(((dFdx3(:,:,2).*dFdx1(:,:,1)) - (dFdx1(:,:,2).*dFdx3(:,:,1)))./((dFdx1(:,:,2).*dFdx2(:,:,1)) - (dFdx2(:,:,2).*dFdx1(:,:,1))));
+tanalpha2=(((dFdx3(:,:,4).*dFdx2(:,:,3)) - (dFdx2(:,:,4).*dFdx3(:,:,3)))./((dFdx2(:,:,4).*dFdx1(:,:,3)) - (dFdx1(:,:,4).*dFdx2(:,:,3))));
+tanbeta2=(((dFdx3(:,:,4).*dFdx1(:,:,3)) - (dFdx1(:,:,4).*dFdx3(:,:,3)))./((dFdx1(:,:,4).*dFdx2(:,:,3)) - (dFdx2(:,:,4).*dFdx1(:,:,3))));
+
+
+end
+
+function [Un_alpha1,Un_alpha2,Un_beta1,Un_beta2]=Uncertainty_propagation_through_mapping_function_gradient(calmat,xg,yg,zg,Uncalcoeff,unwx,unwy,unwz)
+% This function calculates the uncertainty in the mapping function gradients and subsequesntly the angle uncertainty
+%written by Sayantan Bhattacharya January 2016
+[r,c]=size(xg);
+dFdx1=zeros(r,c,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+dFdx2=zeros(r,c,4);
+dFdx3=zeros(r,c,4);
+Un_dFdx1=zeros(r,c,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+Un_dFdx2=zeros(r,c,4);
+Un_dFdx3=zeros(r,c,4);
+
+for gg=1:4
+ a=calmat(:,gg);
+ % Uc=uncoeff(:,gg);
+
+ % Mapping function gradient
+ dFdx1(:,:,gg) = a(2) + 2*a(5)*xg + a(6)*yg + a(8)*zg + 3*a(11)*xg.^2 + 2*a(12)*xg.*yg + ...
+ a(13)*yg.^2 + 2*a(15)*xg.*zg + a(16)*yg.*zg + a(18)*zg.^2;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*xg + 2*a(7)*yg + a(9)*zg + a(12)*xg.^2 + 2*a(13)*xg.*yg + ...
+ 3*a(14)*yg.^2 + a(16)*xg.*zg + 2*a(17)*yg.*zg + a(19)*zg.^2;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*xg + a(9)*yg + 2*a(10)*zg + a(15)*xg.^2 + a(16)*xg.*yg + ...
+ a(17)*yg.^2 + 2*a(18)*xg.*zg + 2*a(19)*yg.*zg;
+ % Second order derivative
+ cfxx=2*a(5) +6*a(11)*xg + 2*a(12)*yg + 2*a(15)*zg ;
+ cfxy=a(6) + 2*a(12)*xg + 2*a(13)*yg + a(16)*zg ;
+ cfxz=a(8)+ 2*a(15)*xg + a(16)*yg + 2*a(18)*zg;
+
+ cfyx=a(6) + 2*a(12)*xg + 2*a(13)*yg + a(16)*zg ;
+ cfyy=2*a(7) +2*a(13)*xg + 6*a(14)*yg + 2*a(17)*zg;
+ cfyz= a(9)+ a(16)*xg + 2*a(17)*yg + 2*a(19)*zg;
+
+ cfzx=a(8)+ 2*a(15)*xg + a(16)*yg + 2*a(18)*zg;
+ cfzy=a(9)+ a(16)*xg + 2*a(17)*yg + 2*a(19)*zg;
+ cfzz=2*a(10)+ 2*a(18)*xg + 2*a(19)*yg;
+
+ %Uncertainty due to Ux, Uy and Uz
+ T1dFdx1=sqrt((cfxx.^2).*(unwx.^2)+(cfxy.^2).*(unwy.^2)+(cfxz.^2).*(unwz.^2));
+ T1dFdx2=sqrt((cfyx.^2).*(unwx.^2)+(cfyy.^2).*(unwy.^2)+(cfyz.^2).*(unwz.^2));
+ T1dFdx3=sqrt((cfzx.^2).*(unwx.^2)+(cfzy.^2).*(unwy.^2)+(cfzz.^2).*(unwz.^2));
+
+
+ % Uncertainty due to cal coefficients
+
+ a=Uncalcoeff(:,gg); % Uncertainty in mapping function coefficients
+
+ T2dFdx1=sqrt(a(2).^2 + (2*a(5)).^2.*xg.^2 + a(6).^2*yg.^2 + a(8).^2*zg.^2 +...
+ (3*a(11)).^2*xg.^4 + (2*a(12)).^2*xg.^2.*yg.^2 + a(13).^2*yg.^4 + ...
+ (2*a(15)).^2*xg.^2.*zg.^2 + a(16).^2*yg.^2.*zg.^2 + a(18).^2*zg.^4);
+
+ T2dFdx2=sqrt(a(3).^2 + a(6).^2*xg.^2 + (2*a(7)).^2*yg.^2 + a(9).^2*zg.^2 +...
+ a(12).^2*xg.^4 + (2*a(13)).^2*xg.^2.*yg.^2 + (3*a(14)).^2*yg.^4 + ...
+ a(16).^2*xg.^2.*zg.^2 + (2*a(17)).^2*yg.^2.*zg.^2 + a(19).^2*zg.^4);
+
+ T2dFdx3=sqrt(a(4).^2 + a(8).^2*xg.^2 + a(9).^2*yg.^2 + (2*a(10)).^2*zg.^2 +...
+ a(15).^2*xg.^4 + a(16).^2*xg.^2.*yg.^2 + a(17).^2*yg.^4 + ...
+ (2*a(18)).^2*xg.^2.*zg.^2 + (2*a(19)).^2*yg.^2.*zg.^2);
+
+ %Total Uncertainty in mapping function gradient
+ Un_dFdx1(:,:,gg)=sqrt(T1dFdx1.^2+T2dFdx1.^2);
+ Un_dFdx2(:,:,gg)=sqrt(T1dFdx2.^2+T2dFdx2.^2);
+ Un_dFdx3(:,:,gg)=sqrt(T1dFdx3.^2+T2dFdx3.^2);
+
+ % Uncertainty in mapping function gradient due to Ux,Uy,Uz
+% Un_dFdx1(:,:,gg)=sqrt(T1dFdx1.^2);
+% Un_dFdx2(:,:,gg)=sqrt(T1dFdx2.^2);
+% Un_dFdx3(:,:,gg)=sqrt(T1dFdx3.^2);
+
+ % Uncertainty in mapping function gradient due to Uai's
+% Un_dFdx1(:,:,gg)=sqrt(T2dFdx1.^2);
+% Un_dFdx2(:,:,gg)=sqrt(T2dFdx2.^2);
+% Un_dFdx3(:,:,gg)=sqrt(T2dFdx3.^2);
+
+end
+% keyboard;
+
+%ANGLE CALCULATION
+N1=((dFdx3(:,:,2).*dFdx2(:,:,1))-(dFdx2(:,:,2).*dFdx3(:,:,1)));
+D1=(dFdx2(:,:,2).*dFdx1(:,:,1)-dFdx1(:,:,2).*dFdx2(:,:,1));
+alpha1 = N1./D1;
+%A1=atan2(N1,D1);
+
+N2=((dFdx3(:,:,4).*dFdx2(:,:,3))-(dFdx2(:,:,4).*dFdx3(:,:,3)));
+D2=(dFdx2(:,:,4).*dFdx1(:,:,3)-dFdx1(:,:,4).*dFdx2(:,:,3));
+alpha2 = N2./D2;
+%A2=atan2(N2,D2);
+
+N3=((dFdx3(:,:,2).*dFdx1(:,:,1))-(dFdx1(:,:,2).*dFdx3(:,:,1)));
+D3=(dFdx1(:,:,2).*dFdx2(:,:,1)-dFdx2(:,:,2).*dFdx1(:,:,1));
+beta1 = N3./D3;
+%B1=atan2(N3,D3);
+
+N4=((dFdx3(:,:,4).*dFdx1(:,:,3))-(dFdx1(:,:,4).*dFdx3(:,:,3)));
+D4=(dFdx1(:,:,4).*dFdx2(:,:,3)-dFdx2(:,:,4).*dFdx1(:,:,3));
+beta2 = N4./D4;
+%B2=atan2(N4,D4);
+
+
+%UNCERTAINTY IN ANGLES
+%uncertainty for alpha1
+mul1=(1./(1+alpha1.^2)).^2;
+term1=(((N1./(D1.^2)).*dFdx2(:,:,2)).^2).*(Un_dFdx1(:,:,1)).^2;
+term2=((dFdx3(:,:,2)./D1 + (N1./D1.^2).*dFdx1(:,:,2)).^2).*(Un_dFdx2(:,:,1)).^2;
+term3=((dFdx2(:,:,2)./D1).^2).*(Un_dFdx3(:,:,1)).^2;
+
+term4=(((N1./(D1.^2)).*dFdx2(:,:,1)).^2).*(Un_dFdx1(:,:,2)).^2;
+term5=((dFdx3(:,:,1)./D1 + (N1./D1.^2).*dFdx1(:,:,1)).^2).*(Un_dFdx2(:,:,2)).^2;
+term6=((dFdx2(:,:,1)./D1).^2).*(Un_dFdx3(:,:,2)).^2;
+
+Un_alpha1=(mul1.*(term1+term2+term3+term4+term5+term6)).^0.5;
+
+%uncertainty for alpha2
+mul2=(1./(1+alpha2.^2)).^2;
+term11=(((N2./(D2.^2)).*dFdx2(:,:,4)).^2).*(Un_dFdx1(:,:,3)).^2;
+term12=((dFdx3(:,:,4)./D2 + (N2./D2.^2).*dFdx1(:,:,4)).^2).*(Un_dFdx2(:,:,3)).^2;
+term13=((dFdx2(:,:,4)./D2).^2).*(Un_dFdx3(:,:,3)).^2;
+
+term14=(((N2./(D2.^2)).*dFdx2(:,:,3)).^2).*(Un_dFdx1(:,:,4)).^2;
+term15=((dFdx3(:,:,3)./D2 + (N2./D2.^2).*dFdx1(:,:,3)).^2).*(Un_dFdx2(:,:,4)).^2;
+term16=((dFdx2(:,:,3)./D2).^2).*(Un_dFdx3(:,:,4)).^2;
+
+Un_alpha2=(mul2.*(term11+term12+term13+term14+term15+term16)).^0.5;
+
+%uncertainty for beta1
+mul3=(1./(1+beta1.^2)).^2;
+term21=(((N3./(D3.^2)).*dFdx1(:,:,2)).^2).*(Un_dFdx2(:,:,1)).^2;
+term22=((dFdx3(:,:,2)./D3 + (N3./D3.^2).*dFdx2(:,:,2)).^2).*(Un_dFdx1(:,:,1)).^2;
+term23=((dFdx1(:,:,2)./D3).^2).*(Un_dFdx3(:,:,1)).^2;
+
+term24=(((N3./(D3.^2)).*dFdx1(:,:,1)).^2).*(Un_dFdx2(:,:,2)).^2;
+term25=((dFdx3(:,:,1)./D3 + (N3./D3.^2).*dFdx2(:,:,1)).^2).*(Un_dFdx1(:,:,2)).^2;
+term26=((dFdx1(:,:,1)./D3).^2).*(Un_dFdx3(:,:,2)).^2;
+
+Un_beta1=(mul3.*(term21+term22+term23+term24+term25+term26)).^0.5;
+
+%uncertainty for beta2
+mul4=(1./(1+beta2.^2)).^2;
+term131=(((N4./(D4.^2)).*dFdx1(:,:,4)).^2).*(Un_dFdx2(:,:,3)).^2;
+term132=((dFdx3(:,:,4)./D4 + (N4./D4.^2).*dFdx2(:,:,4)).^2).*(Un_dFdx1(:,:,3)).^2;
+term133=((dFdx1(:,:,4)./D4).^2).*(Un_dFdx3(:,:,3)).^2;
+
+term134=(((N4./(D4.^2)).*dFdx1(:,:,3)).^2).*(Un_dFdx2(:,:,4)).^2;
+term135=((dFdx3(:,:,3)./D4 + (N4./D4.^2).*dFdx2(:,:,3)).^2).*(Un_dFdx1(:,:,4)).^2;
+term136=((dFdx1(:,:,3)./D4).^2).*(Un_dFdx3(:,:,4)).^2;
+
+Un_beta2=(mul4.*(term131+term132+term133+term134+term135+term136)).^0.5;
+
+end
diff --git a/stereo_uncertainty_codes/stereo_uncertainty_propagation.m b/stereo_uncertainty_codes/stereo_uncertainty_propagation.m
new file mode 100644
index 0000000..08f9bcc
--- /dev/null
+++ b/stereo_uncertainty_codes/stereo_uncertainty_propagation.m
@@ -0,0 +1,290 @@
+function [Un_u,Un_v,Un_w,JU,JV,JW]=stereo_uncertainty_propagation(Unu1,Unv1,Unu2,Unv2,U1,V1,U2,V2,W,Un_alpha1,Un_alpha2,Un_beta1,Un_beta2,tanalpha1,tanalpha2,tanbeta1,tanbeta2,mx,my)
+% This function calculates the uncertainty in the stereo velocity
+% components
+
+
+% Input Variables
+% Nx X Ny grid point array on which individual camera velocity fields are
+% evaluated
+% Unu1= Nx X Ny array of uncertainty in camera 1 U component of velocity using any planar uncertainty method like IM or CS
+% Unv1= Nx X Ny array of uncertainty in camera 1 V component of velocity using any planar uncertainty method like IM or CS
+% Unu2= Nx X Ny array of uncertainty in camera 2 U component of velocity using any planar uncertainty method like IM or CS
+% Unv2= Nx X Ny array of uncertainty in camera 2 V component of velocity using any planar uncertainty method like IM or CS
+% W= Out of plane velocity velocity component
+
+% tanalpha1,tanalpha2,tanbeta1,tanbeta2= tangent of the stereo angles alpha
+% (angle in x-z plane) and beta (angle in y-z plane) for camera 1 and 2
+
+% Un_alpha1,Un_alpha2,Un_beta1,Un_beta2= the uncertainty in stereo
+% angles alpha1, beta1, alpha2, beta2
+
+%The angles and its uncertainties are calculated previously in
+%stereo_angle_uncertainty.m
+
+% mx,my = magnifications in mm(or physical unit)/pixel for the dewarped
+% common grid in x and y direction respectively.
+
+% Output Variables
+% Un_u,Un_v and Un_w are the uncertainties in the stereo velocity
+% components u, v and w.
+
+%JU, JV and JW = Sensitivity coefficients of U, V and W uncertainty propagation equation as given in Table 1 of the stereo uncertainty paper
+
+%written by Sayantan Bhattacharya on 01/12/2016
+
+
+%Get the stereo angles from the tangent of the angles
+A1=atan(tanalpha1);
+A2=atan(tanalpha2);
+B1=atan(tanbeta1);
+B2=atan(tanbeta2);
+
+
+if (max(abs(tanalpha1(:)))+ max(abs(tanalpha2(:))))> (max(abs(tanbeta1(:)))+ max(abs(tanbeta2(:)))) % For horizontal configuration (X-Z plane)
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%Uncertainty horizontal camera%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ %% calculating coefficients for U component
+ %Sensitivity coeeficients for U-component uncertainty propagation
+ Udu1= tanalpha2./(tanalpha2-tanalpha1);
+ Udu2= -tanalpha1./(tanalpha2-tanalpha1);
+ Uda1= ((U1-U2).*sin(A2).*cos(A2))./(sin(A2-A1)).^2;
+ Uda2= ((U2-U1).*sin(A1).*cos(A1))./(sin(A2-A1)).^2;
+
+ lm=zeros(size(Unu1(:)));
+ SigU=zeros(4,4,size(lm,1));
+ JU=zeros(1,4,size(lm,1));
+
+ Un_u=zeros(size(lm));Un_v=zeros(size(lm));Un_w=zeros(size(lm));
+
+ %Calculating U jacobian
+ JU(1,1,:)=Udu1(:); JU(1,2,:)=Udu2(:); JU(1,3,:)=Uda1(:); JU(1,4,:)=Uda2(:);
+
+ %variances or diagonal terms of covariance matrix
+ su1=Unu1(:).^2;
+ su2=Unu2(:).^2;
+ su3=Un_alpha1(:).^2;
+ su4=Un_alpha2(:).^2;
+
+ %correlation coefficients for elemental uncertainties
+ %This is set to zero and thus the covariance terms are neglected for
+ %now
+ ru12=0;ru13=0;ru14=0;ru23=0;ru24=0;ru34=0;
+
+ %Covariance matrix for U
+ SigU(1,1,:)=su1; SigU(1,2,:)=ru12.*sqrt(su1.*su2); SigU(1,3,:)=ru13.*sqrt(su1.*su3); SigU(1,4,:)=ru14.*sqrt(su1.*su4);
+ SigU(2,1,:)=ru12.*sqrt(su1.*su2); SigU(2,2,:)=su2; SigU(2,3,:)=ru23.*sqrt(su2.*su3); SigU(2,4,:)=ru24.*sqrt(su2.*su4);
+ SigU(3,1,:)=ru13.*sqrt(su1.*su3); SigU(3,2,:)=ru23.*sqrt(su2.*su3); SigU(3,3,:)=su3; SigU(3,4,:)=ru34.*sqrt(su3.*su4);
+ SigU(4,1,:)=ru14.*sqrt(su1.*su4); SigU(4,2,:)=ru24.*sqrt(su2.*su4); SigU(4,3,:)=ru34.*sqrt(su3.*su4); SigU(4,4,:)=su4;
+
+
+ %% calculating coefficients for W component
+ %Sensitivity coeeficients for W-component uncertainty propagation
+ Wdu1= 1./(tanalpha2-tanalpha1);
+ Wdu2= -1./(tanalpha2-tanalpha1);
+ Wda1= ((U1-U2).*cos(A2).^2)./(sin(A2-A1)).^2;
+ Wda2= -((U1-U2).*cos(A1).^2)./(sin(A2-A1)).^2;
+
+ SigW=zeros(4,4,size(lm,1));
+ JW=zeros(1,4,size(lm,1));
+
+ %Calculating W jacobian
+ JW(1,1,:)=Wdu1(:); JW(1,2,:)=Wdu2(:); JW(1,3,:)=Wda1(:); JW(1,4,:)=Wda2(:);
+
+ %%variances or diagonal terms of covariance matrix
+ sw1=Unu1(:).^2;
+ sw2=Unu2(:).^2;
+ sw3=Un_alpha1(:).^2;
+ sw4=Un_alpha2(:).^2;
+
+ %correlation coefficients for elemental uncertainties
+ %This is set to zero and thus the covariance terms are neglected for
+ %now
+ rw12=0;rw13=0;rw14=0;rw23=0;rw24=0;rw34=0;
+
+ %Covariance matrix for W
+ SigW(1,1,:)=sw1; SigW(1,2,:)=rw12.*sqrt(sw1.*sw2); SigW(1,3,:)=rw13.*sqrt(sw1.*sw3); SigW(1,4,:)=rw14.*sqrt(sw1.*sw4);
+ SigW(2,1,:)=rw12.*sqrt(sw1.*sw2); SigW(2,2,:)=sw2; SigW(2,3,:)=rw23.*sqrt(sw2.*sw3); SigW(2,4,:)=rw24.*sqrt(sw2.*sw4);
+ SigW(3,1,:)=rw13.*sqrt(sw1.*sw3); SigW(3,2,:)=rw23.*sqrt(sw2.*sw3); SigW(3,3,:)=sw3; SigW(3,4,:)=rw34.*sqrt(sw3.*sw4);
+ SigW(4,1,:)=rw14.*sqrt(sw1.*sw4); SigW(4,2,:)=rw24.*sqrt(sw2.*sw4); SigW(4,3,:)=rw34.*sqrt(sw3.*sw4); SigW(4,4,:)=sw4;
+
+ for i=1:length(lm)
+ Un_u(i)=sqrt(JU(:,:,i)*SigU(:,:,i)*JU(:,:,i)');
+ Un_w(i)=sqrt(JW(:,:,i)*SigW(:,:,i)*JW(:,:,i)');
+ end
+% keyboard;
+ %% calculating coefficients for V component
+ %Sensitivity coeeficients for V-component uncertainty propagation
+ Vdv1= 0.5;
+ Vdv2= 0.5;
+ Vdb1= (mx/my)*(W./2).*(1./(cos(B1).^2));
+ Vdb2= (mx/my)*(W./2).*(1./(cos(B2).^2));
+ VdW = (mx/my)*0.5.*(tanbeta1+tanbeta2);
+
+
+ SigV=zeros(5,5,size(lm,1));
+ JV=zeros(1,5,size(lm,1));
+
+ %Calculating W jacobian
+ JV(1,1,:)=Vdv1(:); JV(1,2,:)=Vdv2(:); JV(1,3,:)=Vdb1(:); JV(1,4,:)=Vdb2(:); JV(1,5,:)=VdW(:);
+
+ %%variances or diagonal terms of covariance matrix
+ sv1=Unv1(:).^2;
+ sv2=Unv2(:).^2;
+ sv3=Un_beta1(:).^2;
+ sv4=Un_beta2(:).^2;
+ sv5=Un_w(:).^2;
+
+ %correlation coefficients for elemental uncertainties
+ %This is set to zero and thus the covariance terms are neglected for
+ %now
+ rv12=0;rv13=0;rv14=0;rv15=0;rv23=0;rv24=0;rv25=0;rv34=0;rv35=0;rv45=0;
+
+ %Covariance matrix for V
+ SigV(1,1,:)=sv1; SigV(1,2,:)=rv12.*sqrt(sv1.*sv2); SigV(1,3,:)=rv13.*sqrt(sv1.*sv3); SigV(1,4,:)=rv14.*sqrt(sv1.*sv4); SigV(1,5,:)=rv15.*sqrt(sv1.*sv5);
+ SigV(2,1,:)=rv12.*sqrt(sv1.*sv2); SigV(2,2,:)=sv2; SigV(2,3,:)=rv23.*sqrt(sv2.*sv3); SigV(2,4,:)=rv24.*sqrt(sv2.*sv4); SigV(2,5,:)=rv25.*sqrt(sv2.*sv5);
+ SigV(3,1,:)=rv13.*sqrt(sv1.*sv3); SigV(3,2,:)=rv23.*sqrt(sv2.*sv3); SigV(3,3,:)=sv3; SigV(3,4,:)=rv34.*sqrt(sv3.*sv4); SigV(3,5,:)=rv35.*sqrt(sv3.*sv5);
+ SigV(4,1,:)=rv14.*sqrt(sv1.*sv4); SigV(4,2,:)=rv24.*sqrt(sv2.*sv4); SigV(4,3,:)=rv34.*sqrt(sv3.*sv4); SigV(4,4,:)=sv4; SigV(4,5,:)=rv45.*sqrt(sv4.*sv5);
+ SigV(5,1,:)=rv15.*sqrt(sv1.*sv5); SigV(5,2,:)=rv25.*sqrt(sv2.*sv5); SigV(5,3,:)=rv35.*sqrt(sv3.*sv5); SigV(5,4,:)=rv45.*sqrt(sv4.*sv5); SigV(5,5,:)=sv5;
+
+ % The uncertainty propagation equation
+ for i=1:length(lm)
+ Un_v(i)=sqrt(JV(:,:,i)*SigV(:,:,i)*JV(:,:,i)');
+ end
+
+ Un_u=reshape(Un_u,size(Unu1,1),size(Unu1,2));
+ Un_v=reshape(Un_v,size(Unu1,1),size(Unu1,2));
+ Un_w=reshape(Un_w,size(Unu1,1),size(Unu1,2));
+
+
+
+elseif (max(abs(tanalpha1(:)))+ max(abs(tanalpha2(:))))< (max(abs(tanbeta1(:)))+ max(abs(tanbeta2(:)))) % For Vertical Configuration (Y-Z plane)
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%Uncertainty vertical camera%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ %% calculating coefficients for V component
+ %Sensitivity coeeficients for V-component uncertainty propagation
+ Vdv1= tanbeta2./(tanbeta2-tanbeta1);
+ Vdv2= -tanbeta1./(tanbeta2-tanbeta1);
+ Vdb1= ((V1-V2).*sin(B2).*cos(B2))./(sin(B2-B1)).^2;
+ Vdb2= ((V2-V1).*sin(B1).*cos(B1))./(sin(B2-B1)).^2;
+
+ lm=zeros(size(Unv1(:)));
+ SigV=zeros(4,4,size(lm,1));
+ JV=zeros(1,4,size(lm,1));
+
+ Un_u=zeros(size(lm));Un_v=zeros(size(lm));Un_w=zeros(size(lm));
+
+ %Calculating V jacobian
+ JV(1,1,:)=Vdv1(:); JV(1,2,:)=Vdv2(:); JV(1,3,:)=Vdb1(:); JV(1,4,:)=Vdb2(:);
+ %keyboard;
+
+ %variances or diagonal terms of covariance matrix
+ sv1=Unv1(:).^2;
+ sv2=Unv2(:).^2;
+ sv3=Un_beta1(:).^2;
+ sv4=Un_beta2(:).^2;
+
+ %correlation coefficients for elemental uncertainties
+ %This is set to zero and thus the covariance terms are neglected for
+ %now
+ rv12=0;rv13=0;rv14=0;rv23=0;rv24=0;rv34=0;
+
+ %Covariance matrix for V
+ SigV(1,1,:)=sv1; SigV(1,2,:)=rv12.*sqrt(sv1.*sv2); SigV(1,3,:)=rv13.*sqrt(sv1.*sv3); SigV(1,4,:)=rv14.*sqrt(sv1.*sv4);
+ SigV(2,1,:)=rv12.*sqrt(sv1.*sv2); SigV(2,2,:)=sv2; SigV(2,3,:)=rv23.*sqrt(sv2.*sv3); SigV(2,4,:)=rv24.*sqrt(sv2.*sv4);
+ SigV(3,1,:)=rv13.*sqrt(sv1.*sv3); SigV(3,2,:)=rv23.*sqrt(sv2.*sv3); SigV(3,3,:)=sv3; SigV(3,4,:)=rv34.*sqrt(sv3.*sv4);
+ SigV(4,1,:)=rv14.*sqrt(sv1.*sv4); SigV(4,2,:)=rv24.*sqrt(sv2.*sv4); SigV(4,3,:)=rv34.*sqrt(sv3.*sv4); SigV(4,4,:)=sv4;
+
+
+
+ %% calculating coefficients for W component
+ %Sensitivity coeeficients for W-component uncertainty propagation
+ Wdv1= 1./(tanbeta2-tanbeta1);
+ Wdv2= -1./(tanbeta2-tanbeta1);
+ Wdb1= ((V1-V2).*cos(B2).^2)./(sin(B2-B1)).^2;
+ Wdb2= -((V1-V2).*cos(B1).^2)./(sin(B2-B1)).^2;
+
+ SigW=zeros(4,4,size(lm,1));
+ JW=zeros(1,4,size(lm,1));
+
+ %Calculating W jacobian
+ JW(1,1,:)=Wdv1(:); JW(1,2,:)=Wdv2(:); JW(1,3,:)=Wdb1(:); JW(1,4,:)=Wdb2(:);
+
+ %%variances or diagonal terms of covariance matrix
+ sw1=Unv1(:).^2;
+ sw2=Unv2(:).^2;
+ sw3=Un_beta1(:).^2;
+ sw4=Un_beta2(:).^2;
+
+ %correlation coefficients for elemental uncertainties
+ %This is set to zero and thus the covariance terms are neglected for
+ %now
+ rw12=0;rw13=0;rw14=0;rw23=0;rw24=0;rw34=0;
+
+
+ %Covariance matrix for W
+ SigW(1,1,:)=sw1; SigW(1,2,:)=rw12.*sqrt(sw1.*sw2); SigW(1,3,:)=rw13.*sqrt(sw1.*sw3); SigW(1,4,:)=rw14.*sqrt(sw1.*sw4);
+ SigW(2,1,:)=rw12.*sqrt(sw1.*sw2); SigW(2,2,:)=sw2; SigW(2,3,:)=rw23.*sqrt(sw2.*sw3); SigW(2,4,:)=rw24.*sqrt(sw2.*sw4);
+ SigW(3,1,:)=rw13.*sqrt(sw1.*sw3); SigW(3,2,:)=rw23.*sqrt(sw2.*sw3); SigW(3,3,:)=sw3; SigW(3,4,:)=rw34.*sqrt(sw3.*sw4);
+ SigW(4,1,:)=rw14.*sqrt(sw1.*sw4); SigW(4,2,:)=rw24.*sqrt(sw2.*sw4); SigW(4,3,:)=rw34.*sqrt(sw3.*sw4); SigW(4,4,:)=sw4;
+
+ for i=1:length(lm)
+ Un_v(i)=sqrt(JV(:,:,i)*SigV(:,:,i)*JV(:,:,i)');
+ Un_w(i)=sqrt(JW(:,:,i)*SigW(:,:,i)*JW(:,:,i)');
+ end
+ %keyboard;
+ %% calculating coefficients for U component
+ %Sensitivity coeeficients for U-component uncertainty propagation
+ Udu1= 0.5;
+ Udu2= 0.5;
+ Uda1= (my/mx)*(W./2).*(1./(cos(A1).^2));
+ Uda2= (my/mx)*(W./2).*(1./(cos(A2).^2));
+ UdW = (my/mx)*0.5.*(tanalpha1+tanalpha2);
+
+
+ SigU=zeros(5,5,size(lm,1));
+ JU=zeros(1,5,size(lm,1));
+
+ %Calculating U jacobian
+ JU(1,1,:)=Udu1(:); JU(1,2,:)=Udu2(:); JU(1,3,:)=Uda1(:); JU(1,4,:)=Uda2(:); JU(1,5,:)=UdW(:);
+
+ %%variances or diagonal terms of covariance matrix
+ su1=Unu1(:).^2;
+ su2=Unu2(:).^2;
+ su3=Un_alpha1(:).^2;
+ su4=Un_alpha2(:).^2;
+ su5=Un_w(:).^2;
+
+ %correlation coefficients for elemental uncertainties
+ %This is set to zero and thus the covariance terms are neglected for
+ %now
+ ru12=0;ru13=0;ru14=0;ru15=0;ru23=0;ru24=0;ru25=0;ru34=0;ru35=0;ru45=0;
+
+ %Couariance matrix for U
+ SigU(1,1,:)=su1; SigU(1,2,:)=ru12.*sqrt(su1.*su2); SigU(1,3,:)=ru13.*sqrt(su1.*su3); SigU(1,4,:)=ru14.*sqrt(su1.*su4); SigU(1,5,:)=ru15.*sqrt(su1.*su5);
+ SigU(2,1,:)=ru12.*sqrt(su1.*su2); SigU(2,2,:)=su2; SigU(2,3,:)=ru23.*sqrt(su2.*su3); SigU(2,4,:)=ru24.*sqrt(su2.*su4); SigU(2,5,:)=ru25.*sqrt(su2.*su5);
+ SigU(3,1,:)=ru13.*sqrt(su1.*su3); SigU(3,2,:)=ru23.*sqrt(su2.*su3); SigU(3,3,:)=su3; SigU(3,4,:)=ru34.*sqrt(su3.*su4); SigU(3,5,:)=ru35.*sqrt(su3.*su5);
+ SigU(4,1,:)=ru14.*sqrt(su1.*su4); SigU(4,2,:)=ru24.*sqrt(su2.*su4); SigU(4,3,:)=ru34.*sqrt(su3.*su4); SigU(4,4,:)=su4; SigU(4,5,:)=ru45.*sqrt(su4.*su5);
+ SigU(5,1,:)=ru15.*sqrt(su1.*su5); SigU(5,2,:)=ru25.*sqrt(su2.*su5); SigU(5,3,:)=ru35.*sqrt(su3.*su5); SigU(5,4,:)=ru45.*sqrt(su4.*su5); SigU(5,5,:)=su5;
+
+ % The uncertainty propagation equation
+ for i=1:length(lm)
+ Un_u(i)=sqrt(JU(:,:,i)*SigU(:,:,i)*JU(:,:,i)');
+ end
+
+ Un_u=reshape(Un_u,size(Unu1,1),size(Unu1,2));
+ Un_v=reshape(Un_v,size(Unu1,1),size(Unu1,2));
+ Un_w=reshape(Un_w,size(Unu1,1),size(Unu1,2));
+
+
+end
+
+
+end
\ No newline at end of file
diff --git a/stereo_uncertainty_codes/uncertainty_in_xyz_calcoeff.m b/stereo_uncertainty_codes/uncertainty_in_xyz_calcoeff.m
new file mode 100644
index 0000000..139b874
--- /dev/null
+++ b/stereo_uncertainty_codes/uncertainty_in_xyz_calcoeff.m
@@ -0,0 +1,921 @@
+function [xg,yg,Uncal,Ux,Uy,Uz,z1grid,caldatamod]=uncertainty_in_xyz_calcoeff(caldata,dispfield,xgrid1,ygrid1)
+% This function calculates the uncertainty in the world coordinates and the
+% uncertainty in the calibration coefficients using uncertainty propagation
+% through triangulation
+% Inputs
+%caldata= calibration job containing calibration information and mapping
+%functions
+%dispfield= disparity vector field U,V based pixel based coordinate grid
+%X,Y
+%xgrid1,ygrid1= dewarped common grid in physical units
+
+%Outputs
+%xg,yg are coordinate grids in physical units where disparity vector is
+%evaluated
+%Uncal= Uncertainty in calibration coefficients
+%Ux,Uy,Uz= Uncertainty in physical coordiantes x,y,z
+%z1grid= mean shift in the triangulated z plane due to any residual
+%disparity, a measure of bias uncertainty in z
+%caldatamod= modified calibration coefficients taking into account the
+%uncertainty in the coefficients
+
+%written by Sayantan Bhattacharya on January 2016
+
+%% Calibration information
+%world coordinates
+allx1data=caldata.allx1data;
+allx2data=caldata.allx2data;
+%camera image coordinates
+allX1data=caldata.allX1data;
+allX2data=caldata.allX2data;
+
+method=caldata.modeltype;
+optionsls=caldata.optionsls;
+
+%polynomial fitting coefficients
+aXcam1=caldata.aXcam1;
+aYcam1=caldata.aYcam1;
+aXcam2=caldata.aXcam2;
+aYcam2=caldata.aYcam2;
+
+calmat=[aXcam1 aYcam1 aXcam2 aYcam2];
+
+%% DISPARITY FIELD information
+Dux=dispfield.U(:,:,1);
+Duy=dispfield.V(:,:,1);%X and Y disparity matrices need to change correlation points referenced to pixel corners to be referenced to pixel centers
+X3=dispfield.X + 0.5;
+Y3=dispfield.Y + 0.5;% correlation X and Y grid points (in pixels)
+Unwx=dispfield.Unwx;% Random uncertainty using image matching in pixels
+Unwy=dispfield.Unwy;
+
+%calculating the length of the common grid in the object plane on which the
+%images are dewarped
+xmin=min(min(xgrid1));
+xmax=max(max(xgrid1));
+ymin=min(min(ygrid1));
+ymax=max(max(ygrid1));
+
+[Imax1,Jmax1]=size(xgrid1);
+% Scaling factor = object plane grid length in x or y(in physical units)mm/ pixel length of images
+scalex=(xmax-xmin)/(Jmax1-1);
+scaley=(ymax-ymin)/(Imax1-1);
+% [scalex scaley]
+%Do coordinate transform between vectors location reported as image indices
+% (From 1 to Imax1) to world coordinates:
+xg = scalex * (X3-1) + xmin; %(in physical units or mm)
+yg = scaley * (Y3-1) + ymin;
+zgrid=zeros(size(xg));
+%% For Bias uncertainty z location
+%scaling the disparity to physical dimensions
+Duxb=Dux.*(scalex);
+Duyb=Duy.*(scaley);
+
+% Removing outliers at the edges
+stdx=rms(Duxb(:)-mean(Duxb(:)));
+stdy=rms(Duyb(:)-mean(Duyb(:)));
+meanx=mean(Duxb(:));
+meany=mean(Duyb(:));
+
+for i=1:size(xg,1)
+ for j=1:size(xg,2)
+
+ if abs(Duxb(i,j))=abs(newzx(lv))
+% % zp(lv)=zbias(lv)+zrand(lv);
+% % zn(lv)=zbias(lv)-zrand(lv);
+% % elseif abs(zg1(lv))=th %&& (max(max(dz)))<=0.1
+
+% Now calculate viewing angles for both linear z and quadratic z
+% calibration models
+
+% Viewing angles for camera 1 calculated at gridpoints xgrid,ygrid;
+
+ aall=[aXcam1 aYcam1 aXcam2 aYcam2];
+ dFdx1=zeros(rows,cols,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+ dFdx2=zeros(rows,cols,4);
+ dFdx3=zeros(rows,cols,4);
+
+ if caldata.modeltype==1
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Mapping the camera coord. to the World Coord. using 1sr order z
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for gg=1:4
+ a=aall(:,gg);
+ dFdx1(:,:,gg) = a(2) + 2*a(5)*xgrid + a(6)*ygrid + a(8)*zgrid + 3*a(10)*xgrid.^2 + ...
+ 2*a(11)*xgrid.*ygrid + a(12)*ygrid.^2 + 2*a(14)*xgrid.*zgrid + a(15)*ygrid.*zgrid;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*xgrid + 2*a(7)*ygrid + a(9)*zgrid + a(11)*xgrid.^2 + ...
+ 2*a(12)*xgrid.*ygrid + 3*a(13)*ygrid.^2 + a(15)*xgrid.*zgrid + 2*a(16)*ygrid.*zgrid;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*xgrid + a(9)*ygrid + a(14)*xgrid.^2 + a(15)*xgrid.*ygrid + a(16)*ygrid.^2;
+ end
+
+ elseif caldata.modeltype==2
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Mapping the camera coord. to the World Coord. using 2nd order z
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for gg=1:4
+ a=aall(:,gg);
+ dFdx1(:,:,gg) = a(2) + 2*a(5).*xgrid + a(6)*ygrid + a(8)*zgrid + 3*a(11)*xgrid.^2 + 2*a(12)*xgrid.*ygrid + ...
+ a(13)*ygrid.^2 + 2*a(15)*xgrid.*zgrid + a(16)*ygrid.*zgrid + a(18)*zgrid.^2;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*xgrid + 2*a(7)*ygrid + a(9)*zgrid + a(12)*xgrid.^2 + 2*a(13)*xgrid.*ygrid + ...
+ 3*a(14)*ygrid.^2 + a(16)*xgrid.*zgrid + 2*a(17)*ygrid.*zgrid + a(19)*zgrid.^2;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*xgrid + a(9)*ygrid + 2*a(10)*zgrid + a(15)*xgrid.^2 + a(16)*xgrid.*ygrid + ...
+ a(17)*ygrid.^2 + 2*a(18)*xgrid.*zgrid + 2*a(19)*ygrid.*zgrid;
+ end
+
+ end
+
+%calculating the viewing angles using formula (7) and (8) from Giordano and Astarita's paper "Spatial resolution of the Stereo PIV technique" (2009)
+
+alphatan(:,:,1)=(((dFdx3(:,:,2).*dFdx2(:,:,1)) - (dFdx2(:,:,2).*dFdx3(:,:,1)))./((dFdx2(:,:,2).*dFdx1(:,:,1)) - (dFdx1(:,:,2).*dFdx2(:,:,1))));
+betatan(:,:,1)=(((dFdx3(:,:,2).*dFdx1(:,:,1)) - (dFdx1(:,:,2).*dFdx3(:,:,1)))./((dFdx1(:,:,2).*dFdx2(:,:,1)) - (dFdx2(:,:,2).*dFdx1(:,:,1))));
+
+%aall=[aXcam1 aYcam1 aXcam2 aYcam2];
+ dFdx1=zeros(rows,cols,4); % the 3rd dimention corresponds to dFdx1 for (X1,Y1,X2,Y2)
+ dFdx2=zeros(rows,cols,4);
+ dFdx3=zeros(rows,cols,4);
+%Xx1c1=aXcam1(2) + 2*aXcam1(5).*x1
+
+
+%Viewing angles for camera 2 calculated at gridpoints x2grid=xgrid+Dux,y2grid=ygrid+Duy
+ if caldata.modeltype==1
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Mapping the camera coord. to the World Coord. using 1sr order z
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for gg=1:4
+ a=aall(:,gg);
+ dFdx1(:,:,gg) = a(2) + 2*a(5)*x2grid + a(6)*y2grid + a(8)*zgrid + 3*a(10)*x2grid.^2 + ...
+ 2*a(11)*x2grid.*y2grid + a(12)*y2grid.^2 + 2*a(14)*x2grid.*zgrid + a(15)*y2grid.*zgrid;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*x2grid + 2*a(7)*y2grid + a(9)*zgrid + a(11)*x2grid.^2 + ...
+ 2*a(12)*x2grid.*y2grid + 3*a(13)*y2grid.^2 + a(15)*x2grid.*zgrid + 2*a(16)*y2grid.*zgrid;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*x2grid + a(9)*y2grid + a(14)*x2grid.^2 + a(15)*x2grid.*y2grid + a(16)*y2grid.^2;
+ end
+
+ elseif caldata.modeltype==2
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Mapping the camera coord. to the World Coord. using 2nd order z
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ for gg=1:4
+ a=aall(:,gg);
+ dFdx1(:,:,gg) = a(2) + 2*a(5).*x2grid + a(6)*y2grid + a(8)*zgrid + 3*a(11)*x2grid.^2 + 2*a(12)*x2grid.*y2grid + ...
+ a(13)*y2grid.^2 + 2*a(15)*x2grid.*zgrid + a(16)*y2grid.*zgrid + a(18)*zgrid.^2;
+
+ dFdx2(:,:,gg) = a(3) + a(6)*x2grid + 2*a(7)*y2grid + a(9)*zgrid + a(12)*x2grid.^2 + 2*a(13)*x2grid.*y2grid + ...
+ 3*a(14)*y2grid.^2 + a(16)*x2grid.*zgrid + 2*a(17)*y2grid.*zgrid + a(19)*zgrid.^2;
+
+ dFdx3(:,:,gg) = a(4) + a(8)*x2grid + a(9)*y2grid + 2*a(10)*zgrid + a(15)*x2grid.^2 + a(16)*x2grid.*y2grid + ...
+ a(17)*y2grid.^2 + 2*a(18)*x2grid.*zgrid + 2*a(19)*y2grid.*zgrid;
+ end
+
+ end
+
+alphatan(:,:,2)=(((dFdx3(:,:,4).*dFdx2(:,:,3)) - (dFdx2(:,:,4).*dFdx3(:,:,3)))./((dFdx2(:,:,4).*dFdx1(:,:,3)) - (dFdx1(:,:,4).*dFdx2(:,:,3))));
+betatan(:,:,2)=(((dFdx3(:,:,4).*dFdx1(:,:,3)) - (dFdx1(:,:,4).*dFdx3(:,:,3)))./((dFdx1(:,:,4).*dFdx2(:,:,3)) - (dFdx2(:,:,4).*dFdx1(:,:,3))));
+%keyboard;
+
+%Display camera angles for reference
+
+% figure(100); subplot(2,2,1);
+% imagesc(atand(alphatan(:,:,1))); colorbar; %caxis([25 30]);
+% title('Camera 1 Angle \alpha1','FontSize',16)
+% subplot(2,2,2);
+% imagesc(atand(alphatan(:,:,2))); colorbar; %caxis([-30 -25]);
+% title('Camera 2 Angle \alpha2','FontSize',16)
+% subplot(2,2,3);
+% imagesc(atand(betatan(:,:,1))); colorbar; %caxis([-2 2]);(end-(zed+10):end-(zed+5))
+% title('Camera 1 Angle \beta1','FontSize',16)
+% subplot(2,2,4);
+% imagesc(atand(betatan(:,:,2))); colorbar; %caxis([-2 2]);
+% title('Camera 1 Angle \beta1','FontSize',16)
+%
+% keyboard;
+
+alpha1=atand(alphatan(:,:,1));
+alpha2=atand(alphatan(:,:,2));
+beta1=atand(betatan(:,:,1));
+beta2=atand(betatan(:,:,2));
+
+% Triangulation using formula (1) of that paper
+
+if max(abs(alpha1(:))+abs(alpha2(:))) > max(abs(beta1(:))+abs(beta2(:)))
+ % if cameras are placed along x axis when alpha>beta
+ if mean(alpha1(:))mean(alpha2(:))
+ %2
+ dz1=(Dux)./(abs(alphatan(:,:,1))+abs(alphatan(:,:,2)));
+ end
+
+elseif max(abs(alpha1(:))+abs(alpha2(:))) < max(abs(beta1(:))+abs(beta2(:)))
+ % if cameras are placed along y axis when alphamean(beta2(:))
+ %4
+ dz1=(Duy)./(abs(betatan(:,:,1))+abs(betatan(:,:,2)));
+ end
+
+end
+
+
+zgrid=zgrid+dz1;% the new z grid
+
+z1grid=zgrid;
+%figure(3);surf(z1grid);title('change in z direction');xlabel('x(mm)');ylabel('y(mm)');
+%figure(4);surf(z1grid);
+
+end
+
+function [a_cam1, a_cam2, aXcam1, aYcam1, aXcam2, aYcam2, convergemessage,Uncal] = fitmodels_mod(allx1data,...
+ allx2data,allX1data,allX2data,modeltype,~)
+% function [a_cam1 a_cam2 aXcam1 aYcam1 aXcam2 aYcam2 convergemessage]=fitcameramodels(allx1data,...
+% allx2data,allX1data,allX2data,modeltype,optionsls)
+% This function ....
+% Inputs:
+% allx1data: Location of target markers (x,y,z) in world coord.
+% for camera 1.
+% allx2data: Location of target markers (x,y,z) in world coord.
+% for camera 2.
+% allX1data: Location of target markers (x,y) in pixel coord.
+% for camera 1.
+% allX2data: Location of target markers (x,y) in pixel coord.
+% for camera 2.
+% Modeltype: A switch to choose which camera model that will be
+% used for the 3D fit. 1 = Cubic XY, Linear Z, 2 =
+% Cubic XY, Quadratic Z, & 3 = Camera Pinhole Model
+% (DLT).
+% optionsls: These are the options for the least squares solver.
+% This variable is determined by optimset
+%
+% Outputs:
+% a_cam1: Polinomial coeff for camera 1 using the pinhole
+% camera model.
+% a_cam2: Polinomial coeff for camera 2 using the pinhole
+% camera model.
+% a_Xcam1: Polinomial coeff for the X component of camera 1
+% using the cubic mapping function.
+% a_Ycam1: Polinomial coeff for the Y component of camera 1
+% using the cubic mapping function.
+% a_Xcam2: Polinomial coeff for the X component of camera 2
+% using the cubic mapping function.
+% a_Ycam2: Polinomial coeff for the Y component of camera 2
+% using the cubic mapping function.
+% convergemessage: Provides a message about the convergence of the
+% non-linear least squares solver. This variable
+% has 4 parts, (X and Y for both cameras). The tag
+% following the Yes or No is from the lsqnonlin
+% function built into matlab. The next to parts
+% are RMS in pixels and the R value.
+
+% Writen by M. Brady
+% Edited and Commented by S. Raben
+
+% This file is part of prana, an open-source GUI-driven program for
+% calculating velocity fields using PIV or PTV.
+%
+% Copyright (C) 2014 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 .
+
+
+% save('allx1dataFLS.mat','allx1data');
+% save('allx2dataFLS.mat','allx2data');
+% save('allX1dataFLS.mat','allX1data');
+% save('allX2dataFLS.mat','allX2data');
+
+% save('allx1datanew.mat','allx1data');
+% save('allx2datanew.mat','allx2data');
+% save('allX1datanew.mat','allX1data');
+% save('allX2datanew.mat','allX2data');
+
+% save('allx1datagood.mat','allx1data');
+% save('allx2datagood.mat','allx2data');
+% save('allX1datagood.mat','allX1data');
+% save('allX2datagood.mat','allX2data');
+
+%keyboard;
+% allx1data(:,1)=(allx1data(:,1)-mean(allx1data(:,1)))./(1*max(abs(allx1data(:,1))));
+% allx1data(:,2)=(allx1data(:,2)-mean(allx1data(:,2)))./(1*max(abs(allx1data(:,2))));
+% allx1data(:,3)=(allx1data(:,3)-mean(allx1data(:,3)))./(1*max(abs(allx1data(:,3))));
+%
+% allx2data(:,1)=(allx2data(:,1)-mean(allx2data(:,1)))./(1*max(abs(allx2data(:,1))));
+% allx2data(:,2)=(allx2data(:,2)-mean(allx2data(:,2)))./(1*max(abs(allx2data(:,2))));
+% allx2data(:,3)=(allx2data(:,3)-mean(allx2data(:,3)))./(1*max(abs(allx2data(:,3))));
+%
+% allX1data(:,1)=(allX1data(:,1)-mean(allX1data(:,1)))./(1*max(abs(allX1data(:,1))));
+% allX1data(:,2)=(allX1data(:,2)-mean(allX1data(:,2)))./(1*max(abs(allX1data(:,2))));
+%
+% allX2data(:,1)=(allX2data(:,1)-mean(allX2data(:,1)))./(1*max(abs(allX2data(:,1))));
+% allX2data(:,2)=(allX2data(:,2)-mean(allX2data(:,2)))./(1*max(abs(allX2data(:,2))));
+
+
+
+
+optionsls=optimset('MaxIter',30000,'MaxFunEvals',30000,'TolX',1e-11,'TolFun',1e-7,...
+ 'LargeScale','off','Display','off','Algorithm','levenberg-marquardt');
+
+
+% Matlab perfers 'Levenberg-Marquardt' over 'levenbergmarquardt' but the
+% optimset function hasn't been updated yet. This should be fixed in later
+% versions of matlab
+% fprintf('LevenbergMarquart warning turned off in fitcameramodel.m\n')
+% warning('off','optim:lsqncommon:AlgorithmConflict')
+
+% Variable Initialization
+a_cam1=[];
+a_cam2=[];
+aXcam1=[];
+aXcam2=[];
+aYcam1=[];
+aYcam2=[];
+
+% This is a switch that will allow for compulation of the camera pinhole
+% coeff. (12) at the same time as the other coeff.. Otherwise it split it
+% into 3 steps.
+% 1 -> for simulatanious compulation
+% 0 -> for split
+allsametime=0;
+
+% Initalization of coeff. for the different models.
+if modeltype==1
+ a=[100; ones(15,1)]; % initial guesss for solver
+elseif modeltype==2
+ a=[100; ones(18,1)]; % initial guesss for solver
+else % camera pinhole model
+
+end
+
+% Defines the percision of the outputs in the fit table in the GUI.
+printformat='%5.4f'; % for converge box message
+printformat1='%7.6f';
+
+% lsqnonlin outputs
+% x,resnorm,residual,exitflag,output,lambda,jacobian
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% If either of the cubic XY camera models are going to be used.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if ((modeltype==1) || (modeltype==2))
+
+ % Set model type
+ alldata.orderz=modeltype;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % get coefficients for X, cam1
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ alldata.allxdata=allx1data;
+ alldata.allXdata=allX1data(:,1);
+ % fit X for cam. 1 to all x data for cam1
+ [aXcam1,resnorm,f_resid,exitflag,~,~,jac1]=lsqnonlin(@(a) poly_3xy_123z(a,alldata),...
+ a,[],[],optionsls); %#ok
+ rmsX1=sqrt(mean((f_resid).^2));%sqrt(resnorm/length(alldata.allXdata));
+ corrcoef=1-sqrt(mean((f_resid./alldata.allXdata).^2));
+ if ((exitflag ~= -1) && (exitflag ~= -2))
+ msgX1=[' yes (' num2str(exitflag) ') ' num2str(rmsX1,printformat) ' ' num2str(corrcoef,printformat1)];
+ else
+ msgX1=[' no (' num2str(exitflag) ') ' num2str(rmsX1,printformat) ' ' num2str(corrcoef,printformat1)];
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % get coefficients for Y, cam1
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ alldata.allxdata=allx1data;
+ alldata.allXdata=allX1data(:,2);
+ % fit Y for cam. 1 to all x data for cam1
+ [aYcam1,resnorm,f_resid,exitflag,~,~,jac2]=lsqnonlin(@(a) poly_3xy_123z(a,alldata),...
+ a,[],[],optionsls); %#ok
+ rmsY1=sqrt(mean((f_resid).^2));%sqrt(resnorm/length(alldata.allXdata));
+ corrcoef=1-sqrt(mean((f_resid./alldata.allXdata).^2));
+ if ((exitflag ~= -1) && (exitflag ~= -2))
+ msgY1=[' yes (' num2str(exitflag) ') ' num2str(rmsY1,printformat) ' ' num2str(corrcoef,printformat1)];
+ else
+ msgY1=[' no (' num2str(exitflag) ') ' num2str(rmsY1,printformat) ' ' num2str(corrcoef,printformat1)];
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % get coefficients for X, cam2
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ alldata.allxdata=allx2data;
+ alldata.allXdata=allX2data(:,1);
+ % fit X for cam. 2 to all x data for cam1
+ [aXcam2,resnorm,f_resid,exitflag,~,~,jac3]=lsqnonlin(@(a) poly_3xy_123z(a,alldata),...
+ a,[],[],optionsls); %#ok
+ rmsX2=sqrt(mean((f_resid).^2));%sqrt(resnorm/length(alldata.allXdata));
+ corrcoef=1-sqrt(mean((f_resid./alldata.allXdata).^2));
+ if ((exitflag ~= -1) && (exitflag ~= -2))
+ msgX2=[' yes (' num2str(exitflag) ') ' num2str(rmsX2,printformat) ' ' num2str(corrcoef,printformat1)];
+ else
+ msgX2=[' no (' num2str(exitflag) ') ' num2str(rmsX2,printformat) ' ' num2str(corrcoef,printformat1)];
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % get coefficients for Y, cam2
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ alldata.allxdata=allx2data;
+ alldata.allXdata=allX2data(:,2);
+ % fit Y for cam. 2 to all x data for cam1
+ [aYcam2,resnorm,f_resid,exitflag,~,~,jac4]=lsqnonlin(@(a) poly_3xy_123z(a,alldata),...
+ a,[],[],optionsls); %#ok
+ rmsY2=sqrt(mean((f_resid).^2));%sqrt(resnorm/length(alldata.allXdata));
+ corrcoef=1-sqrt(mean((f_resid./alldata.allXdata).^2));
+ if ((exitflag ~= -1) && (exitflag ~= -2))
+ msgY2=[' yes (' num2str(exitflag) ') ' num2str(rmsY2,printformat) ' ' num2str(corrcoef,printformat1)];
+ else
+ msgY2=[' no (' num2str(exitflag) ') ' num2str(rmsY2,printformat) ' ' num2str(corrcoef,printformat1)];
+ end
+
+% keyboard;
+
+else
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % If modeltype =3 (camera pinhole, direct linear transformation)
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if allsametime % this way calculates a1 thru a12 all at once
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Determination of Polynomial Coef. for camera 1 using the Pinhole
+ % camera model
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ a=500*ones(11,1);
+ alldata.allxdata=allx1data;
+ alldata.allXdata=allX1data;
+ [a_cam1,resnorm,f_resid,exitflag]=lsqnonlin(@(a) campinhole_LSfun(a,alldata),...
+ a,[],[],optionsls); % fit X for cam. 1 to all x data for cam1
+ %[rows cols]=size(alldata.allxdata);
+ %rmsX=sqrt(resnorm/rows/2);
+ %corrcoef=1-sqrt(mean((f_resid./(reshape(alldata.allXdata,rows*2,1))).^2));
+ %corrcoef=1-sqrt(mean((f_resid./([reshape(alldata.allXdata,rows*2,1)])).^2));
+ Xd=alldata.allXdata(:,1);
+ Yd=alldata.allXdata(:,2);
+ xd=alldata.allxdata(:,1);
+ yd=alldata.allxdata(:,2);
+ zd=alldata.allxdata(:,3);
+ ad=a_cam1;
+ res1=(ad(1).*xd+ad(2).*yd+ad(3).*zd+ad(4))./(ad(9).*xd+ad(10).*yd+ad(11).*zd+1)-Xd;
+ res2=(ad(5).*xd+ad(6).*yd+ad(7).*zd+ad(8))./(ad(9).*xd+ad(10).*yd+ad(11).*zd+1)-Yd;
+ rmsX=sqrt(mean([res1;res2].^2));
+ corrcoef=1-sqrt(mean([res1./Xd;res2./Yd].^2));
+ if ((exitflag ~= -1) && (exitflag ~= -2))
+ msgX1=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ msgY1=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ else
+ msgX1=[' no (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ msgY1=[' no (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ end
+ a_cam1=[a_cam1;1]; % add in a34 (or a12). this was constrained to be 1 in the solver
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Determination of Polynomial Coef. for camera 2 using the Pinhole
+ % camera model
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% a=500*ones(11,1);
+ alldata.allxdata=allx2data;
+ alldata.allXdata=allX2data;
+ [a_cam2,resnorm,f_resid,exitflag]=lsqnonlin(@(a) campinhole_LSfun(a,alldata),...
+ a,[],[],optionsls); % fit X for cam. 1 to all x data for cam1
+ %[rows cols]=size(alldata.allxdata);
+ %rmsX=sqrt(resnorm/rows/2);
+ %corrcoef=1-sqrt(mean((f_resid./(reshape(alldata.allXdata,rows*2,1))).^2));
+ %corrcoef=1-sqrt(mean((f_resid./([reshape(alldata.allXdata,rows*2,1)])).^2));
+ Xd=alldata.allXdata(:,1);
+ Yd=alldata.allXdata(:,2);
+ xd=alldata.allxdata(:,1);
+ yd=alldata.allxdata(:,2);
+ zd=alldata.allxdata(:,3);
+ ad=a_cam2;
+ res1=(ad(1).*xd+ad(2).*yd+ad(3).*zd+ad(4))./(ad(9).*xd+ad(10).*yd+ad(11).*zd+1)-Xd;
+ res2=(ad(5).*xd+ad(6).*yd+ad(7).*zd+ad(8))./(ad(9).*xd+ad(10).*yd+ad(11).*zd+1)-Yd;
+ rmsX=sqrt(mean([res1;res2].^2));
+ corrcoef=1-sqrt(mean([res1./Xd;res2./Yd].^2));
+ if ((exitflag ~= -1) && (exitflag ~= -2))
+ msgX2=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ msgY2=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ else
+ msgX2=[' no (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ msgY2=[' no (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ end
+ a_cam2=[a_cam2;1];
+
+
+ else
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % This way calculates a1 thru a12 on three separate times, they
+ % turn out to be the same
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ a=ones(1,4);
+ alldata.allxdata=allx1data;
+ alldata.allXdata=allX1data(:,1);
+ % fit X for cam. 1 to all x data for cam1
+ [a1_cam1,resnorm1,f_resid1,~]=lsqnonlin(@(a) campinholemod_LSfun(a,alldata),...
+ a,[],[],optionsls);
+
+ a=ones(1,4);
+ alldata.allxdata=allx1data;
+ alldata.allXdata=allX1data(:,2);
+ % fit X for cam. 1 to all y data for cam1
+ [a2_cam1,resnorm2,f_resid2,~]=lsqnonlin(@(a) campinholemod_LSfun(a,alldata),...
+ a,[],[],optionsls);
+
+ a=ones(1,4);
+ alldata.allxdata=allx1data;
+ [rows,~]=size(alldata.allxdata);
+ alldata.allXdata=ones(rows,1);
+ % fit X for cam. 1 to all z data for cam1
+ [a3_cam1,resnorm3,f_resid3,exitflag]=lsqnonlin(@(a) campinholemod_LSfun(a,alldata),...
+ a,[],[],optionsls);
+
+ rmsX=sqrt((resnorm1+resnorm2+resnorm3)/3/rows);
+ corrcoef=1-sqrt(mean([(f_resid1./(allX1data(:,1))).^2;(f_resid2./(allX1data(:,2))).^2;(f_resid3./(ones(rows,1))).^2]));
+
+ msgX1=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ msgY1=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+
+ %[f_resid1 f_resid2 f_resid3]
+
+ a=ones(1,4);
+ alldata.allxdata=allx2data;
+ alldata.allXdata=allX2data(:,1);
+ % fit X for cam. 2 to all x data for cam1
+ [a1_cam2,resnorm1,f_resid1,~]=lsqnonlin(@(a) campinholemod_LSfun(a,alldata),...
+ a,[],[],optionsls);
+
+ a=ones(1,4);
+ alldata.allxdata=allx2data;
+ alldata.allXdata=allX2data(:,2);
+ % fit X for cam. 2 to all y data for cam1
+ [a2_cam2,resnorm2,f_resid2,~]=lsqnonlin(@(a) campinholemod_LSfun(a,alldata),...
+ a,[],[],optionsls);
+
+ a=ones(1,4);
+ alldata.allxdata=allx2data;
+ [rows,~]=size(alldata.allxdata);
+ alldata.allXdata=ones(rows,1);
+ % fit X for cam. 2 to all z data for cam1
+ [a3_cam2,resnorm3,f_resid3,exitflag]=lsqnonlin(@(a) campinholemod_LSfun(a,alldata),...
+ a,[],[],optionsls);
+
+ rmsX=sqrt((resnorm1+resnorm2+resnorm3)/3/rows);
+ corrcoef=1-sqrt(mean([(f_resid1./(allX2data(:,1))).^2;(f_resid2./(allX2data(:,2))).^2;(f_resid3./(ones(rows,1))).^2]));
+
+ msgX2=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+ msgY2=[' yes (' num2str(exitflag) ') ' num2str(rmsX,printformat) ' ' num2str(corrcoef,printformat1)];
+
+ % [f_resid1 f_resid2 f_resid3]
+
+ % Combine the camera coef for each camera
+ a_cam1=[a1_cam1; a2_cam1; a3_cam1];
+ a_cam2=[a1_cam2; a2_cam2; a3_cam2];
+ end
+end
+convergemessage={msgX1;msgY1;msgX2;msgY2};
+% Approximate camera angles calculated for each calibration Mapping
+% function using the ratio of the coefficients.
+% if modeltype==1 || modeltype==2
+%
+% [aXcam1 aYcam1 aXcam2 aYcam2]
+% fprintf('\n Approximate Camera Angles: \n Alpha1 Alpha2 Beta1 Beta2\n');
+% alpha1=atand((aYcam1(4)*aXcam1(3) - aYcam1(3)*aXcam1(4))/(aYcam1(3)*aXcam1(2) - aYcam1(2)*aXcam1(3)));
+% alpha2=atand((aYcam2(4)*aXcam2(3) - aYcam2(3)*aXcam2(4))/(aYcam2(3)*aXcam2(2) - aYcam2(2)*aXcam2(3)));
+% beta1=atand((aYcam1(4)*aXcam1(2) - aYcam1(2)*aXcam1(4))/(aYcam1(2)*aXcam1(3) - aYcam1(3)*aXcam1(2)));
+% beta2=atand((aYcam2(4)*aXcam2(2) - aYcam2(2)*aXcam2(4))/(aYcam2(2)*aXcam2(3) - aYcam2(3)*aXcam2(2)));
+% [alpha1 alpha2 beta1 beta2]
+%
+% elseif modeltype==3
+% [a_cam1;a_cam2]
+% aXcam1=[a_cam1(1,4) a_cam1(1,1) a_cam1(1,2) a_cam1(1,3)]';
+% aYcam1=[a_cam1(2,4) a_cam1(2,1) a_cam1(2,2) a_cam1(2,3)]';
+% aXcam2=[a_cam2(1,4) a_cam2(1,1) a_cam2(1,2) a_cam2(1,3)]';
+% aYcam2=[a_cam2(2,4) a_cam2(2,1) a_cam2(2,2) a_cam2(2,3)]';
+% [aXcam1 aYcam1 aXcam2 aYcam2]
+%
+% fprintf('\n Approximate Camera Angles: \n Alpha1 Alpha2 Beta1 Beta2\n');
+% alpha1=atand((aYcam1(4)*aXcam1(3) - aYcam1(3)*aXcam1(4))/(aYcam1(3)*aXcam1(2) - aYcam1(2)*aXcam1(3)));
+% alpha2=atand((aYcam2(4)*aXcam2(3) - aYcam2(3)*aXcam2(4))/(aYcam2(3)*aXcam2(2) - aYcam2(2)*aXcam2(3)));
+% beta1=atand((aYcam1(4)*aXcam1(2) - aYcam1(2)*aXcam1(4))/(aYcam1(2)*aXcam1(3) - aYcam1(3)*aXcam1(2)));
+% beta2=atand((aYcam2(4)*aXcam2(2) - aYcam2(2)*aXcam2(4))/(aYcam2(2)*aXcam2(3) - aYcam2(3)*aXcam2(2)));
+% [alpha1 alpha2 beta1 beta2]
+%
+% aXcam1=[];aYcam1=[];aXcam2=[];aYcam2=[];
+%
+% end
+% save('aXcam1.mat','aXcam1');
+% save('aXcam2.mat','aXcam2');
+% save('aYcam1.mat','aYcam1');
+% save('aYcam2.mat','aYcam2');
+% Turning the warning back on
+%warning('on','optim:lsqncommon:AlgorithmConflict')
+
+% Calculating fit coefficient Uncertainty
+Np=size(allx1data,1);
+Na=length(aXcam1);
+Ndof=Np-Na;
+
+SigresX1=sqrt(Np/Ndof)*rmsX1*eye(Na);
+SigresY1=sqrt(Np/Ndof)*rmsY1*eye(Na);
+SigresX2=sqrt(Np/Ndof)*rmsX2*eye(Na);
+SigresY2=sqrt(Np/Ndof)*rmsY2*eye(Na);
+
+% Cov1=inv(jac1'*jac1);
+% Cov2=inv(jac2'*jac2);
+% Cov3=inv(jac3'*jac3);
+% Cov4=inv(jac4'*jac4);
+
+Cov1=(jac1'*jac1);
+Cov2=(jac2'*jac2);
+Cov3=(jac3'*jac3);
+Cov4=(jac4'*jac4);
+
+Uxcam1=sqrt(diag(Cov1\SigresX1.^2));
+Uycam1=sqrt(diag(Cov2\SigresY1.^2));
+Uxcam2=sqrt(diag(Cov3\SigresX2.^2));
+Uycam2=sqrt(diag(Cov4\SigresY2.^2));
+
+Uncal=[Uxcam1 Uycam1 Uxcam2 Uycam2];
+
+end