pivchal05b_ixx_postprocess.m

% uncertainty IXX testing for 2005B case

clear all;
addpath Z:\Planar_Uncertainty_work\codes\export_fig ;
istring1=sprintf(['%%s%%0%0.0fd.','mat'],3);
istring2=sprintf(['%%s%%0%0.0fd.','txt'],3);%sprintf('%%s%%0%0.0fd',3);
savematdir='Z:\Planar_Uncertainty_work\Results\Matfiles\06_18_2015\';


% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\2005B\';
% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\2005B\newixx\';
outdir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\2005B\ws32IM\';

biasfact=0;

if ~exist(outdir1,'dir')
    mkdir(outdir1);
end
casename1='cclsq2_';
outdir=[outdir1,casename1]; 
dstore=0;
printfig=0;

% dirscc='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\2005B\PIV_scc_chal05_diaCClsq_pass4_';
% dirscc='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\2005B\newtest\PIV_scc_chal05_diaCClsq_pass4_';
% dirscc2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\2005B\PIV_scc_chal05_diaCClsq_pass5_';

% dirscc='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2005B\PIV_scc_chal05_diaCClsq_pass4_';

% dirscc='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2005B\newixx\PIV_scc_chal05_diaCClsq_pass4_';
% dirscc2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2005B\PIV_scc_chal05_diaCClsq_pass5_';
% dirscc='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\2005B\PIV_scc_chal05_diaCClsq_pass4_';
% dirscc2='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\2005B\PIV_scc_chal05_diaCClsq_pass5_';

% dirscc='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\2005B\PIV_scc_chal05_diaCClsq_pass4_';
% dirscc2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\2005B\ws32IM\PIV_scc_chal05_diaCClsq_pass4_';
dirscc='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\2005B\PIV_scc_chal05_diaCClsq_pass4_';
dirscc2='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\2005B\PIV_scc_chal05_diaCClsq_pass4_';


% TRUE SOLUTION DIRECTORY
dirtrue='Z:\My_planar_uncertainty_work_and_other_stuff_that_I_need_to_sort\Uncertainty Work\2005B\resultsixx\newresults\CaseBt';

f=[9 29 49 69 89 109];
D=[2.8 2.8];
toff=0;
for k=1:6
    
    n=f(k+toff);

sccsol=load(sprintf(istring1,dirscc,n));
sccsol2=load(sprintf(istring1,dirscc2,n));

truesol=load(sprintf(istring2,dirtrue,n+1));


%Getting True Solution
x=truesol(:,1);
y=truesol(:,2);
u=truesol(:,3);
v=truesol(:,4);

Xt=zeros(41,88);
Yt=zeros(41,88);
Ut=zeros(41,88);
Vt=zeros(41,88);
Nx=41;Ny=88;
m=1;
for i=1:41
    for j=1:88
        Xt(i,j)=x(m);
        Yt(i,j)=y(m);
        Ut(i,j)=u(m);
        Vt(i,j)=v(m);
        m=m+1;
    end
end

Utrue(:,:,k)=Ut;
Vtrue(:,:,k)=Vt;

% keyboard;
%Getting Solution

Us(:,:,k)=(sccsol.U(:,:,1)./2);
Vs(:,:,k)=(sccsol.V(:,:,1)./2);
Eval(:,:,k)=sccsol.Eval(:,:,1);


MI(:,:,k)=reshape(sccsol.uncertainty(:,5),Nx,Ny);
IXX2(:,:,k)=reshape(sccsol2.uncertainty(:,9)./2,Nx,Ny);
IYY2(:,:,k)=reshape(sccsol2.uncertainty(:,10)./2,Nx,Ny);
% IXX2(:,:,k)=reshape(sccsol3.uncertainty(:,25)./2,Nx,Ny);
% IYY2(:,:,k)=reshape(sccsol3.uncertainty(:,26)./2,Nx,Ny);

Autod(:,:,k)=reshape(sccsol.uncertainty(:,6),Nx,Ny);

Ixx=reshape((sccsol.uncertainty(:,7))./2,Nx,Ny);
Iyy=reshape((sccsol.uncertainty(:,8))./2,Nx,Ny);

Gpx(:,:,k)=reshape(sccsol.uncertainty(:,15),Nx,Ny);
Gpy(:,:,k)=reshape(sccsol.uncertainty(:,16),Nx,Ny);
Spx(:,:,k)=biasfact*reshape(sccsol.uncertainty(:,17),Nx,Ny);
Spy(:,:,k)=biasfact*reshape(sccsol.uncertainty(:,18),Nx,Ny);

%  %correct gradient (if desired)
[Udiff]=socdiff(Us(:,:,k),16,1);
[Vdiff]=socdiff(Vs(:,:,k),16,2);
% Ixxt = sqrt(Ixx.^2 - (D(1)^2/16)*(Udiff).^2 );
% Iyyt = sqrt(Iyy.^2 - (D(2)^2/16)*(Vdiff).^2 );
% Ixxt=Ixx;
% Iyyt=Iyy;
Ixxt= real(sqrt(Ixx.^2 - (Autod(:,:,k).^2/16).*(Udiff).^2 ));
Iyyt= real(sqrt(Iyy.^2 - (Autod(:,:,k).^2/16).*(Vdiff).^2 ));

%Bias error
mewx(:,:,k)=(abs(Spx(:,:,k)-Gpx(:,:,k)))./2;
mewy(:,:,k)=(abs(Spy(:,:,k)-Gpy(:,:,k)))./2;

%Scaling for ixx
Ixxt=Ixxt./sqrt(MI(:,:,k));
Iyyt=Iyyt./sqrt(MI(:,:,k));
% Final second order moment uncertainty
IXX1(:,:,k)=sqrt(Ixxt.^2+mewx(:,:,k).^2);
IYY1(:,:,k)=sqrt(Iyyt.^2+mewy(:,:,k).^2);
% IXX1(:,:,k)=sqrt(Ixxt.^2);
% IYY1(:,:,k)=sqrt(Iyyt.^2);

err_sccu(:,:,k)=abs(Utrue(:,:,k)-Us(:,:,k));
err_sccv(:,:,k)=abs(Vtrue(:,:,k)-Vs(:,:,k));

noabserrx(:,:,k)=(Utrue(:,:,k)-Us(:,:,k));
noabserry(:,:,k)=(Vtrue(:,:,k)-Vs(:,:,k));

th=0.2;
tmp_eu=err_sccu(:,:,k);
tmp_ev=err_sccv(:,:,k);
tmp_eu(tmp_eu<th)=0;
tmp_ev(tmp_ev<th)=0;
Evalu(:,:,k)=tmp_eu;
Evalv(:,:,k)=tmp_ev;

% % rms_dx(k)=rms(reshape(err_sccu(:,:,k),Nx*Ny,1));
% rms_dy(k)=rms(reshape(err_sccv(:,:,k),Nx*Ny,1));
% % rms_ixx(k)=rms(reshape(IXX1(:,:,k),Nx*Ny,1));
% rms_iyy(k)=rms(reshape(IYY1(:,:,k),Nx*Ny,1));
% 
% tex=err_sccu(:,:,k);
% tixx=IXX1(:,:,k);
% tex=tex(:);
% tixx=tixx(:);
% 
% % iex=find(tex>1);
% % tex(iex)=[];
% % tixx(iex)=[];
% 
% rms_dx(k)=rms(tex);
% rms_ixx(k)=rms(tixx);


% k1=k1+1;
% clear Us Vs Ut Vt MI Ixx Iyy Ixxd Iyyd;
% % figure;imagesc(Us);
% % figure;imagesc(Ut);
% % % figure;imagesc(Vs);
% % % figure;imagesc(Vt);
% % % %figure;hist(MI(:),100);
% keyboard;

% x_centroid(:,k1)=sccsol.uncertainty(:,21);
% y_centroid(:,k1)=sccsol.uncertainty(:,22);
% x_centroid_scc(:,k1)=sccsol.uncertainty(:,23);
% y_centroid_scc(:,k1)=sccsol.uncertainty(:,24);
if dstore==1
% Saving Diameters
DCCx(:,:,k)=reshape(sccsol.uncertainty(:,21),Nx,Ny);
DCCy(:,:,k)=reshape(sccsol.uncertainty(:,22),Nx,Ny);
DGccx(:,:,k)=reshape(sccsol.uncertainty(:,29),Nx,Ny);
DGccy(:,:,k)=reshape(sccsol.uncertainty(:,30),Nx,Ny);
Dauto1x(:,:,k)=reshape(sccsol.uncertainty(:,23),Nx,Ny);
Dauto1y(:,:,k)=reshape(sccsol.uncertainty(:,24),Nx,Ny);
Dauto2x(:,:,k)=reshape(sccsol.uncertainty(:,27),Nx,Ny);
Dauto2y(:,:,k)=reshape(sccsol.uncertainty(:,28),Nx,Ny);
D1x(:,:,k)=reshape(sccsol.uncertainty(:,31),Nx,Ny);
D1y(:,:,k)=reshape(sccsol.uncertainty(:,32),Nx,Ny);
end
% keyboard;
end
% keyboard;
% save('2005Brmserror_prana_valid72.mat','rms_dx','rms_dy','rms_ixx','rms_iyy');
% keyboard;

% figure;hist(MI(:),100);
% dxz=find(IXX3(:)==-10);
% dyz=find(IYY3(:)==-10);
%[length(dxz) length(dyz)]

% clim=0.04;
% figure;
% % subplot(2,3,1);imagesc(err_sccu(:,:,1));caxis([0 clim]);title('error at t=9');
% % subplot(2,3,2);imagesc(err_sccu(:,:,3));caxis([0 clim]);title('error at t=49');
% % subplot(2,3,3);imagesc(err_sccu(:,:,5));caxis([0 clim]);title('error at t=89');
% % subplot(2,3,4);imagesc(IXX1(:,:,1));caxis([0 clim]);title('Unx IXX at t=9');
% % subplot(2,3,5);imagesc(IXX1(:,:,3));caxis([0 clim]);title('Unx IXX at t=49');
% % subplot(2,3,6);imagesc(IXX1(:,:,5));caxis([0 clim]);title('Unx IXX at t=89');
% 
% subplot(2,3,1);imagesc(err_sccu(:,:,1));caxis([0 clim]);title('error at t=10');
% subplot(2,3,2);imagesc(err_sccu(:,:,2));caxis([0 clim]);title('error at t=30');
% subplot(2,3,3);imagesc(err_sccu(:,:,3));caxis([0 clim]);title('error at t=50');
% subplot(2,3,4);imagesc(err_sccu(:,:,4));caxis([0 clim]);title('error at t=70');
% subplot(2,3,5);imagesc(err_sccu(:,:,5));caxis([0 clim]);title('error at t=90');
% subplot(2,3,6);imagesc(err_sccu(:,:,6));caxis([0 clim]);title('error at t=110');
% keyboard;
% % clim=3;
% % figure;
% % subplot(2,3,1);imagesc(Evalu(:,:,1));caxis([0 clim]);title('error at t=10');
% % subplot(2,3,2);imagesc(Evalu(:,:,2));caxis([0 clim]);title('error at t=30');
% % subplot(2,3,3);imagesc(Evalu(:,:,3));caxis([0 clim]);title('error at t=50');
% % subplot(2,3,4);imagesc(Evalu(:,:,4));caxis([0 clim]);title('error at t=70');
% % subplot(2,3,5);imagesc(Evalu(:,:,5));caxis([0 clim]);title('error at t=90');
% % subplot(2,3,6);imagesc(Evalu(:,:,6));caxis([0 clim]);title('error at t=110');
% % 
% % clim=3;
% % figure;
% % subplot(2,3,1);imagesc(Evalv(:,:,1));caxis([0 clim]);title('error at t=10');
% % subplot(2,3,2);imagesc(Evalv(:,:,2));caxis([0 clim]);title('error at t=30');
% % subplot(2,3,3);imagesc(Evalv(:,:,3));caxis([0 clim]);title('error at t=50');
% % subplot(2,3,4);imagesc(Evalv(:,:,4));caxis([0 clim]);title('error at t=70');
% % subplot(2,3,5);imagesc(Evalv(:,:,5));caxis([0 clim]);title('error at t=90');
% % subplot(2,3,6);imagesc(Evalv(:,:,6));caxis([0 clim]);title('error at t=110');

%%
noabserrx=noabserrx(:);
noabserry=noabserry(:);
err_sccu=err_sccu(:);
err_sccv=err_sccv(:);
IXX1=IXX1(:);
IYY1=IYY1(:);
IXX2=IXX2(:);
IYY2=IYY2(:);

% IXX2=abs(IXX2);
% IYY2=abs(IYY2);

MI=MI(:);
if dstore==1
   DCCx=DCCx(:);
   DCCy=DCCy(:);
   DGccx=DGccx(:);
   DGccy=DGccy(:);
   Dauto1x=Dauto1x(:);
   Dauto1y=Dauto1y(:);
   Dauto2x=Dauto2x(:);
   Dauto2y=Dauto2y(:);
   D1x=D1x(:);
   D1y=D1y(:);
end
% DGccx=DGccx./sqrt(MI);
% DGccy=DGccy./sqrt(MI);

% save([savematdir,'2005BSOM.mat'],'err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2');

% save('2005BSOMnew_toterr.mat','err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2');
% keyboard;
%%
%zero value removal
inx=(find(IXX1(:)==0));
iny=(find(IYY1(:)==0));
% keyboard;
IXX1(inx)=[];
IYY1(iny)=[];
IXX2(inx)=[];
IYY2(iny)=[];
err_sccu(inx)=[];
err_sccv(iny)=[];
noabserrx(inx)=[];
noabserry(iny)=[];
% MI(inx)=[];
%%
%nan value removal
inx2=(find(isnan(IXX2)));
iny2=(find(isnan(IXX2)));

%keyboard;
IXX1(inx2)=[];
IYY1(iny2)=[];
IXX2(inx2)=[];
IYY2(iny2)=[];
err_sccu(inx2)=[];
err_sccv(iny2)=[];

inx3=(find(isnan(IXX1)));
iny3=(find(isnan(IXX1)));
IXX1(inx3)=[];
IYY1(iny3)=[];
IXX2(inx3)=[];
IYY2(iny3)=[];
err_sccu(inx3)=[];
err_sccv(iny3)=[];

[length(inx) length(iny) length(inx2) length(iny2) length(inx3) length(iny3)]
%%
%invalid vector removal
elm=0.2;%0.2;
inx1=(find(err_sccu(:)>elm));
iny1=(find(err_sccv(:)>elm));
[length(inx1) length(iny1)]
%keyboard;
IXX1(inx1)=[];
IYY1(iny1)=[];
IXX2(inx1)=[];
IYY2(iny1)=[];
err_sccu(inx1)=[];
err_sccv(iny1)=[];

if dstore==1
    DCCx(inx1)=[];
    DGccx(inx1)=[];
    Dauto1x(inx1)=[];
    Dauto2x(inx1)=[];
    D1x(inx1)=[];
    
    DCCy(iny1)=[];
    DGccy(iny1)=[];
    Dauto1y(iny1)=[];
    Dauto2y(iny1)=[];
    D1y(iny1)=[];
end
noabserrx(inx1)=[];
noabserry(iny1)=[];
%%
figure;hist(noabserrx,30);
figure;hist(noabserry,30);

[mean(noabserrx) mean(noabserry)]

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%keyboard;

gmx=(find(err_sccu(:)<=1));
gmy=(find(err_sccv(:)<=1));
ngmx=length(gmx);
ngmy=length(gmy);

%keyboard;
% N=r*s*t;
% [(N-ngmx)*100/N (N-ngmy)*100/N];

cnt=0;
for j=1:length(IXX1)
if err_sccu(j)<=IXX1(j)
cnt=cnt+1;
end
end
covx=100*cnt/ngmx;

cnt=0;
for j=1:length(IYY1)
if err_sccv(j)<=IYY1(j)
cnt=cnt+1;
end
end
covy=100*cnt/ngmy;

cnt=0;
for j=1:length(IXX2)
if err_sccu(j)<=IXX2(j)
cnt=cnt+1;
end
end
covx2=100*cnt/ngmx;

cnt=0;
for j=1:length(IYY2)
if err_sccv(j)<=IYY2(j)
cnt=cnt+1;
end
end
covy2=100*cnt/ngmy;

 [covx covx2]
 [covy covy2]
% save('2005BSOM_new.mat','err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2','covx','covy','covx2','covy2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
lw=3;
fs=18;
%%
% % vec=linspace(0,0.2,60);
% % 
% % Nrex=histc(err_sccu,vec);
% % Nrey=histc(err_sccv,vec);
% % Nrdx=histc(IXX1,vec);
% % Nrdy=histc(IYY1,vec);
% % Nrdx2=histc(IXX2,vec);
% % Nrdy2=histc(IYY2,vec);
% % 
% % ll=1:100:15e2;%max(Nrdx);
% % ll2=1:100:20e2;%max(Nrdy);
% % 
% % sigu=rms(err_sccu(err_sccu<=0.15)).*ones(length(ll),1);
% % sigv=rms(err_sccv(err_sccv<=0.15)).*ones(length(ll2),1);
% % 
% % mixx1=median(IXX1(IXX1<=0.15)).*ones(length(ll),1);
% % miyy1=median(IYY1(IYY1<=0.15)).*ones(length(ll2),1);
% % 
% % mixx2=median(IXX2(IXX2<=0.15)).*ones(length(ll),1);
% % miyy2=median(IYY2(IYY2<=0.15)).*ones(length(ll2),1);
% % 
% % 
% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
% % plot(vec,Nrex,'b-o',vec,Nrdx,'r-*',vec,Nrdx2,'m-o',sigu,ll,'k--',mixx1,ll,'r--',mixx2,ll,'m--');
% % %legend('|errorx|','IXX AMM','IXX DMM','rms|errorx|');
% % legend('|e_x|','I_X_X','\delta_x','\sigma_{|e_x|}','\sim{I_X_X}','\sim{\delta_x}');
% % title('Error and Uncertainty histogram');
% % %title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)]));
% % ylabel('No. of Data Points');xlabel('pixels');
% % 
% % figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
% % plot(vec,Nrey,'b-o',vec,Nrdy,'r-*',vec,Nrdy2,'m-o',sigv,ll2,'k--',miyy1,ll2,'r--',miyy2,ll2,'m--');
% % %legend('|Errory|','IYY AMM','IYY DMM','rms|errory|');
% % legend('|e_y|','I_Y_Y','\delta_y','\sigma_{|e_y|}','\sim{I_Y_Y}','\sim{\delta_y}');
% % title('Error and Uncertainty histogram');
% % %title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)]));
% % ylabel('No. of Data Points');xlabel('pixels');
%%
thresh=0.1;
vec=linspace(0,thresh,40);
% vec2=linspace(0,0.05,100);
% vec3=linspace(0.03,0.1,100);
% vec4=linspace(0,0.03,100);
% vec5=linspace(0.06,0.1,100);

Nrex=histc(err_sccu,vec);
Nrey=histc(err_sccv,vec);
Nrdx=histc(IXX1,vec);
Nrdy=histc(IYY1,vec);
Nrdx2=histc(IXX2,vec);
Nrdy2=histc(IYY2,vec);

ll=1:100:15e2;%max(Nrdx);
ll2=1:100:20e2;%max(Nrdy);
% ll=1:100:15e1;%max(Nrdx);
% ll2=1:100:20e1;%max(Nrdy);

sigu=rms(err_sccu(err_sccu<=thresh)).*ones(length(ll),1);
sigv=rms(err_sccv(err_sccv<=thresh)).*ones(length(ll2),1);

% mixx1=median(IXX1(IXX1<=0.1)).*ones(length(ll),1);
% miyy1=median(IYY1(IYY1<=0.1)).*ones(length(ll2),1);
% 
% mixx2=median(IXX2(IXX2<=0.1)).*ones(length(ll),1);
% miyy2=median(IYY2(IYY2<=0.1)).*ones(length(ll2),1);

mixx1=rms(IXX1(IXX1<=thresh)).*ones(length(ll),1);
miyy1=rms(IYY1(IYY1<=thresh)).*ones(length(ll2),1);

mixx2=rms(IXX2(IXX2<=thresh)).*ones(length(ll),1);
miyy2=rms(IYY2(IYY2<=thresh)).*ones(length(ll2),1);

fprintf(['RMSerr X','\n']);
fprintf(num2str([sigu(1) mixx1(1) mixx2(1)],'%02.4f\t'));
fprintf(['\nRMSerr Y','\n']);
fprintf(num2str([sigv(1) miyy1(1) miyy2(1)],'%02.4f\t'));
fprintf('\n');


figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
plot(vec,Nrex,'k-',vec,Nrdx,'r-',vec,Nrdx2,'r--',sigu,ll,'k-',mixx1,ll,'r-',mixx2,ll,'r--');

hold off;
%legend('|errorx|','IXX AMM','IXX DMM','rms|errorx|');
%legend('|e_x|','I_X_X','\delta_x','\sigma_{|e_x|}','\sim{I_X_X}','\sim{\delta_x}');
% title('Error and Uncertainty histogram(X)');
ylabel('No. of Data Points');xlabel('pixels');
box on;
title(sprintf(['Ixx=',num2str(covx),'  IMx=',num2str(covx2)]));
if printfig==1
    print(gcf,'-dpng',[outdir,'rmsUhist.png']);
end
% export_fig(gcf,[outdir,'rmsUhist.png']);

figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);hold on;
plot(vec,Nrey,'k-',vec,Nrdy,'b-',vec,Nrdy2,'b--',sigv,ll2,'k-',miyy1,ll2,'b-',miyy2,ll2,'b--');
% plot(vec2,Nppryl,'color',rgb('DarkCyan'));
% plot(vec3,Nppryh,'color',rgb('DarkCyan'));
% plot(vec4,Nmiyl,'color',rgb('SkyBlue'));
% plot(vec5,Nmiyh,'color',rgb('SkyBlue'));
hold off;
%legend('|Errory|','IYY AMM','IYY DMM','rms|errory|');
%legend('|e_y|','I_Y_Y','\delta_y','\sigma_{|e_y|}','\sim{I_Y_Y}','\sim{\delta_y}');
% title('Error and Uncertainty histogram(Y)');
ylabel('No. of Data Points');xlabel('pixels');
box on;
title(sprintf(['Iyy=',num2str(covy),'  IMy=',num2str(covy2)]));
if printfig==1
    print(gcf,'-dpng',[outdir,'rmsVhist.png']);
end
    %%
Nbin=20;

maxixx=0.08;
minixx=0;
vu=linspace(minixx,maxixx,Nbin);
maxiyy=0.08;
miniyy=0;
vv=linspace(miniyy,maxiyy,Nbin);

rmserru1=zeros(1,length(vu));
rmserrv1=zeros(1,length(vv));
getixx1=zeros(1,length(vu));
getiyy1=zeros(1,length(vv));
numcx1=zeros(1,length(vu));
numcy1=zeros(1,length(vv));
covxbin1=zeros(1,length(vu));
covybin1=zeros(1,length(vv));

rmserru2=zeros(1,length(vu));
rmserrv2=zeros(1,length(vv));
getixx2=zeros(1,length(vu));
getiyy2=zeros(1,length(vv));
numcx2=zeros(1,length(vu));
numcy2=zeros(1,length(vv));
covxbin2=zeros(1,length(vu));
covybin2=zeros(1,length(vv));

%p=1;

for p=1:length(vu)
    if p<length(vu)
        [val1]= find((IXX1(:)<vu(p+1)) & (IXX1(:)>=vu(p)));
        [val2]= find((IXX2(:)<vu(p+1)) & (IXX2(:)>=vu(p)));
    elseif p==length(vu)
        [val1]= find((IXX1(:)>=vu(p)));
        [val2]= find((IXX2(:)>=vu(p)));
    end
   numcx1(p)=(length(val1));
   numcx2(p)=(length(val2));
   tmp1=(err_sccu(val1));
   tmp2=(err_sccu(val2));

   tempixx1=IXX1(val1);
   tempixx2=IXX2(val2);
   cvrx1=find(tmp1<tempixx1);
   cvrx2=find(tmp2<tempixx2);
   covxbin1(p)=100*length(cvrx1)/numcx1(p);
   covxbin2(p)=100*length(cvrx2)/numcx2(p);
   rmserru1(p)=rms(tmp1);
   rmserru2(p)=rms(tmp2);
   getixx1(p)=rms(IXX1(val1));
   getixx2(p)=rms(IXX2(val2));
    
end
% getixx(length(vu))=length(IXX1)-sum(numcx);
% rmserru(length(vu))=rms(err_sccu(val));
val1=0;val2=0;tmp1=0;tmp2=0;
%p=1;

%axis([0 0.4 0 0.4]);

for p=1:length(vv)
    if p<length(vv)
        [val1]= find((IYY1(:)<vv(p+1)) & (IYY1(:)>=vv(p)));
        [val2]= find((IYY2(:)<vv(p+1)) & (IYY2(:)>=vv(p)));
    elseif p==length(vv)
        [val1]= find((IYY1(:)>=vv(p)));
        [val2]= find((IYY2(:)>=vv(p)));
    end
   numcy1(p)=(length(val1));
   numcy2(p)=(length(val2));
   tmp1=(err_sccv(val1));
   tmp2=(err_sccv(val2));

   tempiyy1=IYY1(val1);
   tempiyy2=IYY2(val2);
   cvry1=find(tmp1<tempiyy1);
   cvry2=find(tmp2<tempiyy2);
   covybin1(p)=100*length(cvry1)/numcy1(p);
   covybin2(p)=100*length(cvry2)/numcy2(p);
 
   rmserrv1(p)=rms(tmp1);
   rmserrv2(p)=rms(tmp2);
   getiyy1(p)=rms(IYY1(val1));
   getiyy2(p)=rms(IYY2(val2));

    
end


figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
plot(getixx1(1:end-1),rmserru1(1:end-1),'r-');hold on;
plot(getixx2(1:end-1),rmserru2(1:end-1),'r--');
plot(getiyy1(1:end-1),rmserrv1(1:end-1),'b-');
plot(getiyy2(1:end-1),rmserrv2(1:end-1),'b--');
%plot(getixx(1:end-1),getixx(1:end-1),'k--');
%plot(getixx1(1:end-1),getixx1(1:end-1),'k--');
plot(getiyy1(1:end-1),getiyy1(1:end-1),'k--');hold off;
title('Rms error vs Uncertainty');
%legend('Ixx vs rms|errorx|','Iyy vs rms|errory|','Ixx=rms|errorx|');
ylabel('rms_{|e|}');xlabel('uncertainty(pixels)');
if printfig==1
    print(gcf,'-dpng',[outdir,'rmserror.png']);

end
%keyboard;

figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
plot(getixx1(1:end-1),covxbin1(1:end-1),'r-');hold on;
plot(getixx2(1:end-1),covxbin2(1:end-1),'r--');
plot(getiyy1(1:end-1),covybin1(1:end-1),'b-');
plot(getiyy2(1:end-1),covybin2(1:end-1),'b--');
%plot(getixx(1:end-1),getixx(1:end-1),'k--');
plot(getiyy1(1:end-1),68.5*ones(length(vv)-1),'k--');hold off;
title('Coverage in each Uncertainty bin');
%legend('Ixx','Iyy','\sigma=68.5%');
ylabel('Coverage');xlabel('uncertainty(pixels)');
if printfig==1
    print(gcf,'-dpng',[outdir,'coverage.png']);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Diameter histograms

if dstore==1
    
    vecdx=linspace(1,6,60);
    
    Nccx=histc(DCCx,vecdx);
    Ngccx=histc(DGccx,vecdx);
    Nauto1x=histc(Dauto1x,vecdx);
    Nauto2x=histc(Dauto2x,vecdx);
    Nd1x=histc(D1x,vecdx);
    
    vecdy=linspace(1,6,60);
    
    Nccy=histc(DCCy,vecdy);
    Ngccy=histc(DGccy,vecdy);
    Nauto1y=histc(Dauto1y,vecdy);
    Nauto2y=histc(Dauto2y,vecdy);
    Nd1y=histc(D1y,vecdy);
    
    figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
    hold on;
    plot(vecdx,Nccx,'r-o');
    plot(vecdx,Ngccx,'b-o');
    plot(vecdx,Nauto1x,'g-o');
    plot(vecdx,Nauto2x,'m-o');
    plot(vecdx,Nd1x,'c-o');
    hold off;
    xlabel('Diameter (in pixels)');
    ylabel('No. of points');
    title('X Diameter');
    legend('CC','Gcc','Auto1','Auto2','AvgAuto');
%     print(gcf,'-dpng',[outdir,'Xdia.png']);
    
    figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
    hold on;
    plot(vecdy,Nccy,'r-o');
    plot(vecdy,Ngccy,'b-o');
    plot(vecdy,Nauto1y,'g-o');
    plot(vecdy,Nauto2y,'m-o');
    plot(vecdy,Nd1y,'c-o');
    hold off;
    xlabel('Diameter (in pixels)');
    ylabel('No. of points');
    title('Y Diameter');
    legend('CC','Gcc','Auto1','Auto2','AvgAuto');
%     print(gcf,'-dpng',[outdir,'Ydia.png']);
    
    
end
%% MI Histogram

% vecmi=linspace(0,80,100);
% Nmi=histc(MI(:),vecmi);
% figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
% plot(vecmi,Nmi);
% xlabel('MI');
% ylabel('No. of points');
% title('MI histogram');
% print(gcf,'-dpng',[outdir,'MIhist.png']);