vortexring_original_uncertainty_postprocess.m

% Plot vortex ring cam1 planar uncertainty using second order moment

clear all;
 lw=3;fs=18;
 
% dirsol='Z:\Planar_Uncertainty_work\Vortex_ring_cam1_planar_uncertainty\PIV_scc_pass4_';
% dirsol='Z:\Planar_Uncertainty_work\Results\Ixx_with_cclsq_dia\vortex_ring\PIV_scc_pass4_';
% dirsol='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\Vortex_Ring\PIV_scc_pass4_';
% dirsol='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\With_zero_center_gauss\Vortex_Ring\PIV_scc_pass4_';


% dirsol='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\Vortex_Ring\PIV_scc_pass4_';
% % dirsol2='Z:\Planar_Uncertainty_work\Vortex_ring_cam1_planar_uncertainty\PIV_scc_pass5_';
% dirsol2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\Vortex_Ring\PIV_scc_pass5_';

% dirsol='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Vortex_Ring\PIV_scc_pass4_';
% dirsol='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Vortex_Ring\newixx\PIV_scc_pass4_';
% dirsol2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Vortex_Ring\PIV_scc_pass5_';
% dirsol='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\Vortex_Ring\PIV_scc_pass4_';
% dirsol2='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\Vortex_Ring\PIV_scc_pass5_';

% dirsol='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Vortex_Ring\PIV_scc_pass4_';
% dirsol2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Vortex_Ring\ws32IM\PIV_scc_pass4_';

dirsol='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\Vortex_Ring\PIV_scc_pass4_';
dirsol2='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\Vortex_Ring\PIV_scc_pass4_';


dirtrue='Z:\Stereo_Uncertainty_Work\VortexRing_True_solution\MATLAB_Data\Eval_1\E_eval1_TrueSolution_';
savematdir='Z:\Planar_Uncertainty_work\Results\Matfiles\06_18_2015\';
printfig=0;

% outdir1='Z:\Planar_Uncertainty_work\Results\vortex_ring_solution\plots\original\';
% outdir1='Z:\Planar_Uncertainty_work\Results\Ixx_with_cclsq_dia\plots\Vortexring\';
% outdir1='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\plots\Vortex_Ring\';
% outdir1='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\With_zero_center_gauss\plots\Vortex_Ring\thresh1\';
% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\Vortex_Ring\newixx\';
outdir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\Vortex_Ring\ws32IM\';

biasfact=0;

if ~exist(outdir1,'dir')
    mkdir(outdir1);
end
casename1='diacclsq';
outdir=[outdir1,casename1]; 


stereojob=load('Z:\Stereo_Uncertainty_Work\vortex_ring_stereo_uncertainty_analysis\Camera1Camera3\Jobs\StereoJob_cam13_fulljob.mat');
stereosol=load('Z:\Stereo_Uncertainty_Work\vortex_ring_stereo_uncertainty_analysis\Camera1Camera3\Reconstructed_fields\Willert\Camera1Camera3_3Cfields\piv_2d3c_cam1cam3_pass_4_00025.mat');

xscale=1/11.9690;%stereojob.datasave.scaling.wil.xscale;
yscale=1/11.9653;%stereojob.datasave.scaling.wil.yscale;
xorigin=576.8623;
yorigin=395.8564;
soltemp=load([dirsol,num2str(30,'%05.0f'),'.mat']);
Xsol1=soltemp.X;
Ysol1=soltemp.Y;

[Xt,Yt]=meshgrid(-34:0.45:20,-27:0.45:27);

Xt1=Xt./xscale+xorigin;
Yt1=Yt./yscale+yorigin;

ax=min(Xt1(:));
bx=max(Xt1(:));
ay=min(Yt1(:));
by=max(Yt1(:));

xvec=Xsol1(1,:);
yvec=Ysol1(:,1);

xind=find(xvec>ax & xvec<bx);
yind=find(yvec>ay & yvec<by);

Xs=Xsol1(yind,xind);
Ys=Ysol1(yind,xind);

% keyboard;
% % xvec=(linspace(1,1024,127)-576.8623)./11.9690;
% % yvec=(linspace(1,1024,127)-395.8564)./11.9653;
% % 
% % [Xs1,Ys1]=meshgrid(xvec,yvec);
% % yind=10:89;xind=22:101;
% % 
% % Xs=Xs1(yind,xind);
% % Ys=Ys1(yind,xind);

% Xs=1000.*stereosol.X;
% Ys=1000.*stereosol.Y;




spf=1;%0.001;

Sx=127;
Sy=127;


offset=24;
offset1=24;
N=50;

% for i=1:N
%     TS=load([dirtrue,num2str(i+offset,'%05.0f'),'.mat']);
%     Ut(:,:,i)=TS.velocity_data.x_velocity;
%     Vt(:,:,i)=TS.velocity_data.y_velocity;
%     
% end
% keyboard;
for i=1:N
    TS=load([dirtrue,num2str(i+offset,'%05.0f'),'.mat']);
    sol=load([dirsol,num2str(i+offset1,'%05.0f'),'.mat']);
    sol5=load([dirsol2,num2str(i+offset1,'%05.0f'),'.mat']);
    % True Solution
    uttemp=TS.velocity_data.x_velocity;
    vttemp=TS.velocity_data.y_velocity;
    %Interpolate true solution onto solution grid
%     ut=interp2(Xt,Yt,uttemp,Xs,Ys,'cubic',0);
%     vt=interp2(Xt,Yt,vttemp,Xs,Ys,'cubic',0);
%     ut=interp2(Xt1,Yt1,uttemp,Xs,Ys,'cubic',0);
%     vt=interp2(Xt1,Yt1,vttemp,Xs,Ys,'cubic',0);
    ut=interp2(Xt1,Yt1,uttemp,Xs,Ys,'linear',0);
    vt=interp2(Xt1,Yt1,vttemp,Xs,Ys,'linear',0);
    %convert to pix/frame
    Ut(:,:,i)=(spf/xscale).*ut;
    Vt(:,:,i)=(spf/yscale).*vt;
    %planar processing soluiton
%     Us(:,:,i)=sol.U(10:89,22:101);
%     Vs(:,:,i)=sol.V(10:89,22:101);
% keyboard;
    Us(:,:,i)=sol.U(yind,xind);
    Vs(:,:,i)=sol.V(yind,xind);
%     keyboard;
    %error calculation
    error_U(:,:,i)=abs(Ut(:,:,i)-Us(:,:,i));
    error_V(:,:,i)=abs(Vt(:,:,i)-Vs(:,:,i));
    
    noabserrx(:,:,i)=(Ut(:,:,i)-Us(:,:,i));
    noabserry(:,:,i)=(Vt(:,:,i)-Vs(:,:,i));

    
%     %error calculation
%     error_U1(:,:,i)=(Ut(:,:,i)-Us(:,:,i));
%     error_V1(:,:,i)=(Vt(:,:,i)-Vs(:,:,i));
    
    % Image Matching uncertainty
    timx=reshape(sol5.uncertainty(:,9),Sx,Sy);
    Imx(:,:,i)=timx(yind,xind);
    timy=reshape(sol5.uncertainty(:,10),Sx,Sy);
    Imy(:,:,i)=timy(yind,xind);
    tnim=reshape(sol5.uncertainty(:,11),Sx,Sy);
    Nim(:,:,i)=tnim(yind,xind);
    
    % Storing Ixx and other parameters.
    tppr=reshape(sol.uncertainty(:,1),Sx,Sy);
    PPR(:,:,i)=tppr(yind,xind);
    
    tmi=reshape(sol.uncertainty(:,5),Sx,Sy);
    MI(:,:,i)=tmi(yind,xind);
    
    Autod=reshape(sol.uncertainty(:,6),Sx,Sy);
    Autod=Autod(yind,xind);
    Ixx=reshape(sol.uncertainty(:,7),Sx,Sy);
    Ixx=Ixx(yind,xind);
    Iyy=reshape(sol.uncertainty(:,8),Sx,Sy);
    Iyy=Iyy(yind,xind);
    
    %Gcc peak location  
    Gpx=reshape(sol.uncertainty(:,15),Sx,Sy);
    Gpx=Gpx(yind,xind);
    Gpy=reshape(sol.uncertainty(:,16),Sx,Sy);
    Gpy=Gpy(yind,xind);
    %Scc peak location
    Spx=reshape(sol.uncertainty(:,17),Sx,Sy);
    Spx=biasfact*Spx(yind,xind);
    Spy=reshape(sol.uncertainty(:,18),Sx,Sy);
    Spy=biasfact*Spy(yind,xind);
    %10:89,22:101;
    %Bias error
    mewx=abs(Spx-Gpx);
    mewy=abs(Spy-Gpy);
    % correct for gradient
    [Udiff]=socdiff(Us(:,:,i),8,1);
    [Vdiff]=socdiff(Vs(:,:,i),8,2);
   
    Ixx= real(sqrt(Ixx.^2 - (Autod.^2/16).*(Udiff).^2 ));
    Iyy= real(sqrt(Iyy.^2 - (Autod.^2/16).*(Vdiff).^2 ));
    
    %Scaling for ixx
    Ixx=Ixx./sqrt(MI(:,:,i));
    Iyy=Iyy./sqrt(MI(:,:,i));
    
%     Ixxtemp(:,:,i)=Ixx;
%     Iyytemp(:,:,i)=Iyy;
    
    % Final second order moment uncertainty
    IXX1(:,:,i)=sqrt(Ixx.^2+mewx.^2);
    IYY1(:,:,i)=sqrt(Iyy.^2+mewy.^2);
%     IXX1(:,:,i)=sqrt(Ixx.^2);
%     IYY1(:,:,i)=sqrt(Iyy.^2);
    
end
% keyboard;
% %% plot error and uncertainty distribution
% 
% sigx=IXX1(:);%Ixxtemp(:);
% sigy=IYY1(:);%Iyytemp(:);
% 
% xx=linspace(-0.3,0.3,60);
% Ixg=zeros(length(xx),length(sigx));
% 
% for i=1:length(sigx)
%     Ixg(:,i)=squeeze(normpdf(xx,0,sigx(i))');
% end
% IXG=nanmean(Ixg,2);
% % IXG=sum(Ixg,2);
% 
% 
% Nerrx=histc(error_U1(:),xx);
% Nixxg=histc(IXG(:),xx);
% 
% figure;plot(xx,IXG,'b-',xx,Nerrx,'k-');
%%
% % figure;imagesc(Ut(:,:,1));caxis([-1.5 3.5]);%caxis([-0.05 0.15]);
% % figure;imagesc((spf/xscale).*stereosol.U);caxis([-1.5 3.5]);%caxis([-0.05 0.15]);
% % 
% % figure;imagesc(Us(:,:,2));caxis([-1.5 3.5]);
% % 
% % figure;imagesc(Us(:,:,1)-Ut(:,:,1));caxis([-0.5 0.5]);
% % 
% % figure;hist(error_U(:),100);
% % figure;hist(error_V(:),100);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % figure(1);set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
% % subplot(1,2,1);imagesc(Us(:,:,i));caxis([-0.5 2]);
% % title('U solution');
% % xlabel('X');ylabel('Y');axis image tight;colorbar;
% % subplot(1,2,2);imagesc(Vs(:,:,i));caxis([-0.5 2]);
% % title('V solution');
% % xlabel('X');ylabel('Y');axis image tight;colorbar;
% % print(gcf,'-dpng',[outdir1,'Solution.png']);
% % 
% % figure(2);set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
% % subplot(1,2,1);imagesc(Ut(:,:,i));caxis([-0.5 2]);
% % title('U true');xlabel('X');ylabel('Y');axis image tight;colorbar;
% % subplot(1,2,2);imagesc(Vt(:,:,i));caxis([-0.5 2]);
% % title('V true');xlabel('X');ylabel('Y');axis image tight;colorbar;
% % print(gcf,'-dpng',[outdir1,'True.png']);
% % 
% % % for i=1:N
% %     figure(11);set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
% %     subplot(1,2,1);imagesc(mean(error_U,3));caxis([0 0.5]);
% %     title('mean U error');xlabel('X');ylabel('Y');axis image tight;colorbar;
% %     subplot(1,2,2);imagesc(mean(error_V,3));caxis([0 0.5]);
% %     title('mean V error');xlabel('X');ylabel('Y');axis image tight;colorbar;
% %     print(gcf,'-dpng',[outdir1,'error.png']);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     pause(0.2);
% end
% keyboard;
%% convert to vectors

noabserrx=noabserrx(:);
noabserry=noabserry(:);

IXX1=IXX1(:);
IYY1=IYY1(:);
Imx=Imx(:);
Imy=Imy(:);
error_U=error_U(:);
error_V=error_V(:);

IXX2=Imx;
IYY2=Imy;

% save([savematdir,'VRSOM.mat'],'error_U','error_V','IXX1','IYY1','IXX2','IYY2');
% 
% keyboard;
% save('VRSOMnew_toterr.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2');
% keyboard;
%% zero value removal
inx=(find(IXX1(:)==0));
iny=(find(IYY1(:)==0));
%keyboard;
IXX1(inx)=[];
IYY1(iny)=[];
Imx(inx)=[];
Imy(iny)=[];
error_U(inx)=[];
error_V(iny)=[];
noabserrx(inx)=[];
noabserry(iny)=[];
%%
%nan value removal
inx2=(find(isnan(Imx)));
iny2=(find(isnan(Imy)));

%keyboard;
IXX1(inx2)=[];
IYY1(iny2)=[];
Imx(inx2)=[];
Imy(iny2)=[];
error_U(inx2)=[];
error_V(iny2)=[];

inx3=(find(isnan(IXX1)));
iny3=(find(isnan(IXX1)));
IXX1(inx3)=[];
IYY1(iny3)=[];
Imx(inx3)=[];
Imy(iny3)=[];
error_U(inx3)=[];
error_V(iny3)=[];

[length(inx) length(iny) length(inx2) length(iny2) length(inx3) length(iny3)]
%% invalid vector removal
thresh1=1;
inx1=(find(error_U(:)>thresh1));
iny1=(find(error_V(:)>thresh1));
%keyboard;
IXX1(inx1)=[];
IYY1(iny1)=[];
Imx(inx1)=[];
Imy(iny1)=[];
error_U(inx1)=[];
error_V(iny1)=[];
noabserrx(inx1)=[];
noabserry(iny1)=[];
%%
figure;hist(noabserrx,60);
figure;hist(noabserry,60);

[mean(noabserrx) mean(noabserry)]
%% Find Coverage

gmx=(find(error_U(:)<=1));
gmy=(find(error_V(:)<=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 error_U(j)<=IXX1(j)
cnt=cnt+1;
end
end
covx=100*cnt/ngmx;

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

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

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

[covx covx2]
[covy covy2]
 
% [rms(error_U) rms(error_V)]

% IXX2=Imx;
% IYY2=Imy;

% save('VRSOM.mat','error_U','error_V','IXX1','IYY1','IXX2','IYY2','covx','covy','covx2','covy2');

 %% Plot error and uncertainty histogram
thresh=0.2;
 
vec=linspace(0,thresh,60);

Nrex=histc(error_U,vec);
Nrey=histc(error_V,vec);
Nrdx=histc(IXX1,vec);
Nrdy=histc(IYY1,vec);
Nrdx2=histc(Imx,vec);
Nrdy2=histc(Imy,vec);


ll=1:100:max(Nrex);
ll2=1:100:max(Nrey);


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

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

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

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;

% title('Error and Uncertainty histogram(X)');
title(sprintf(['Ixx=',num2str(covx),'  IMx=',num2str(covx2)]));
%title(sprintf(['2003B Ixx(SCC)','coverage=',num2str(covx)]));
ylabel('No. of Data Points');xlabel('pixels');
axis([0 0.2 0 inf]);
if printfig==1
    print(gcf,'-dpng',[outdir,'Uhist2.png']);
end

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--');
hold off;

% title('Error and Uncertainty histogram(Y)');
title(sprintf(['Iyy=',num2str(covy),'  IMy=',num2str(covy2)]));
%title(sprintf(['2003B Iyy(SCC)','coverage=',num2str(covy)]));
ylabel('No. of Data Points');xlabel('pixels');
axis([0 0.2 0 inf]);
if printfig==1
    print(gcf,'-dpng',[outdir,'Vhist2.png']);
end

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');
%%
plotflag=1;

if plotflag==1;
    
    IXX2=Imx;
    IYY2=Imy;
    err_sccu=error_U;
    err_sccv=error_V;
    
    Nbin=15;
    
    maxixx=0.1;
    minixx=0;
    vu=linspace(minixx,maxixx,Nbin);
    maxiyy=0.1;
    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));
        %    tmp1(tmp1>1)=[];
        %    tmp2(tmp2>1)=[];
        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)=sqrt(sum((tmp1-mean(tmp1)).^2)/numcx1(p));
        %    rmserru2(p)=sqrt(sum((tmp2-mean(tmp2)).^2)/numcx2(p));
        %    rmserru1(p)=rms(tmp1-mean(tmp1));
        %    rmserru2(p)=rms(tmp2-mean(tmp2));
        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));
        %    tmp1(tmp1>1)=[];
        %    tmp2(tmp2>1)=[];
        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)=sqrt(sum((tmp1-mean(tmp1)).^2)/numcy1(p));
        %    rmserrv2(p)=sqrt(sum((tmp2-mean(tmp2)).^2)/numcy2(p));
        %    rmserrv1(p)=rms(tmp1-mean(tmp1));
        %    rmserrv2(p)=rms(tmp2-mean(tmp2));
        rmserrv1(p)=rms(tmp1);
        rmserrv2(p)=rms(tmp2);
        getiyy1(p)=rms(IYY1(val1));
        getiyy2(p)=rms(IYY2(val2));
        
        
    end
    % getixx(length(vu))=length(IXX1)-sum(numcx);
    % rmserru(length(vu))=rms(err_sccu(val));
%     figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
%     plot(getixx1(1:end-1),100.*numcx1(1:end-1)./sum(numcx1),'r-');hold on;
%     plot(getixx2(1:end-1),100.*numcx2(1:end-1)./sum(numcx2),'r--');
%     plot(getiyy1(1:end-1),100.*numcy1(1:end-1)./sum(numcy1),'b-');
%     plot(getiyy2(1:end-1),100.*numcy2(1:end-1)./sum(numcy2),'b--');hold off;
%     title('% of points in each Uncertainty bin');
%     %legend('Ixx Histogram','Iyy Histogram');
%     ylabel('% of total points');xlabel('uncertainty(pixels)');
%     print(gcf,'-dpng',[outdir1,'number.png']);
    
    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,'rms error in binned uncertainty.png']);
    end
%     print(gcf,'-dpng',[outdir1,'rms error in binned uncertainty.png']);
    
    %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('Covergae in each Uncertainty bin');
    %legend('Ixx','Iyy','\sigma=68.5%');
    ylabel('Coverage');xlabel('uncertainty(pixels)');
    if printfig==1
        print(gcf,'-dpng',[outdir,'covg in each bin.png']);
    end
%     print(gcf,'-dpng',[outdir1,'covg in each bin.png']);
end