postprocess_prana_2003B.m

%analyzing PIV challenge 2003B first case.
clear all;
biasfact=0;
dstore=0;
printfig=0;
%% Input Directories
% Ixx directroy
% basedir1='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\2003B\';
% basedir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2003B\';
% % basedir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2003B\newixx\';
% % dirixx=[basedir1,'PIV_scc_finalchal03_pass4_'];
% % % IM directroy
% % % basedir2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\2003B\';
% % basedir2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\2003B\';
% % dirim=[basedir2,'PIV_scc_finalchal03_pass5_'];

% basedir1='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\2003B\';
% dirixx=[basedir1,'PIV_scc_finalchal03_pass4_'];
% % IM directroy
% % basedir2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\2003B\';
% basedir2='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\2003B\';
% dirim=[basedir2,'PIV_scc_finalchal03_pass5_'];

basedir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\2003B\';
dirixx=[basedir1,'PIV_scc_finalchal03_pass4_'];
% IM directroy
% basedir2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\2003B\';
basedir2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\2003B\ws32IM\';
dirim=[basedir2,'PIV_scc_finalchal03_pass4_'];


% basedir1='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\2003B\';
% dirixx=[basedir1,'PIV_scc_finalchal03_pass4_'];
% % IM directroy
% % basedir2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\2003B\';
% basedir2='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\2003B\';
% dirim=[basedir2,'PIV_scc_finalchal03_pass4_'];

%% Output Directory

savematdir='Z:\Planar_Uncertainty_work\Results\Matfiles\06_18_2015\';
if ~exist(savematdir,'dir')
    mkdir(savematdir);
end
% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\2003B\';
% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\2003B\newixx\';
outdir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\2003B\ws32IM\';

if ~exist(outdir1,'dir')
    mkdir(outdir1);
end
casename1='cclsq';
outdir=[outdir1,casename1]; 
%% True Solution directory
truesol=load('Z:\Planar_Uncertainty_work\Results\True_solution\2003B\solutiononly.mat');
[r,s,t]=size(truesol.u);
% t=2;
%% Loop through solution and True solution 
for i=1:t
    

sccsol= load([dirixx,num2str(2*i-1,'%03.0f'),'.mat']);
sccsol2=load([dirim,num2str(2*i-1,'%03.0f'),'.mat']);

Us=sccsol.U(end:-1:1,:,1);
Vs=sccsol.V(end:-1:1,:,1);

Ut=truesol.u(:,:,i)';
Vt=truesol.v(:,:,i)';

% figure;imagesc(Us);figure;imagesc(Ut);
% figure;imagesc(Vs);figure;imagesc(Vt);
%keyboard;
x=sccsol.X(:);
y=sccsol.Y(:);
a=sort(unique(x));
b=sort(unique(y));
N=length(x);

for n=1:N
        I=find(b==y(n));
        J=find(a==x(n));
        autod(I,J)=sccsol.uncertainty(n,6);
%         autod(I,J)=mean([sccsol.uncertainty(n,21) sccsol.uncertainty(n,22)]);
        np(I,J)=sccsol2.uncertainty(n,11);
        imx(I,J)=sccsol2.uncertainty(n,9);
        imy(I,J)=sccsol2.uncertainty(n,10);
        gpx(I,J)=sccsol.uncertainty(n,15);
        gpy(I,J)=sccsol.uncertainty(n,16);
        spx(I,J)=sccsol.uncertainty(n,17);
        spy(I,J)=sccsol.uncertainty(n,18);
        
        ppr(I,J)=sccsol.uncertainty(n,1);
        prmsr(I,J)=sccsol.uncertainty(n,2);
        pce(I,J)=sccsol.uncertainty(n,3);
        entropy(I,J)=sccsol.uncertainty(n,4);
        mi(I,J)=sccsol.uncertainty(n,5);
        
        Ixx1(I,J)=sccsol.uncertainty(n,7);
        Iyy1(I,J)=sccsol.uncertainty(n,8);
        %Peak2(I,J)=snr(n,5);
end
Gpx(:,:,i)=gpx(end:-1:1,:);
Gpy(:,:,i)=gpy(end:-1:1,:);
Spx(:,:,i)=biasfact*spx(end:-1:1,:);
Spy(:,:,i)=biasfact*spy(end:-1:1,:);

Ppr(:,:,i)=ppr(end:-1:1,:);
Prmsr(:,:,i)=prmsr(end:-1:1,:);
Pce(:,:,i)=pce(end:-1:1,:);
Entropy(:,:,i)=entropy(end:-1:1,:);
MI(:,:,i)=mi(end:-1:1,:);
Np(:,:,i)=np(end:-1:1,:);


IXXd(:,:,i)=imx(end:-1:1,:);
IYYd(:,:,i)=imy(end:-1:1,:);
Autod(:,:,i)=autod(end:-1:1,:);

err_sccu(:,:,i)=abs(Us-Ut);
err_sccv(:,:,i)=abs(Vs-Vt);

noabserrx(:,:,i)=(Ut-Us);
noabserry(:,:,i)=(Vt-Vs);

Ixx(:,:,i)=Ixx1(end:-1:1,:);
Iyy(:,:,i)=Iyy1(end:-1:1,:);

[Udiff]=socdiff(Us,8,1);
[Vdiff]=socdiff(Vs,8,2);

if min(Ixx(:))<0
keyboard;
end
if min(Iyy(:))<0
keyboard;
end

% correct gradient (if desired)
% Ixx(:,:,i) = real(sqrt(Ixx(:,:,i).^2 - (D(1)^2/16)*(Udiff).^2 ));
% Iyy(:,:,i) = real(sqrt(Iyy(:,:,i).^2 - (D(2)^2/16)*(Vdiff).^2 ));

Ixx(:,:,i) = real(sqrt(Ixx(:,:,i).^2 - (Autod(:,:,i).^2/16).*(Udiff).^2 ));
Iyy(:,:,i) = real(sqrt(Iyy(:,:,i).^2 - (Autod(:,:,i).^2/16).*(Vdiff).^2 ));

% Ixx(:,:,i) = real(sqrt(Ixx(:,:,i).^2 + (Autod(:,:,i).^2/16).*(Udiff).^2 ));
% Iyy(:,:,i) = real(sqrt(Iyy(:,:,i).^2 + (Autod(:,:,i).^2/16).*(Vdiff).^2 ));

% Ixx(:,:,i)=sqrt(Ixx(:,:,i).^2+(Spx(:,:,i)-Gpx(:,:,i)).^2);
% Iyy(:,:,i)=sqrt(Iyy(:,:,i).^2+(Spy(:,:,i)-Gpy(:,:,i)).^2);
Ny=length(a);Nx=length(b);


if dstore==1
% Saving Diameters
DCCx(:,:,i)=reshape(sccsol.uncertainty(:,21),Nx,Ny);
DCCy(:,:,i)=reshape(sccsol.uncertainty(:,22),Nx,Ny);
DGccx(:,:,i)=reshape(sccsol.uncertainty(:,29),Nx,Ny);
DGccy(:,:,i)=reshape(sccsol.uncertainty(:,30),Nx,Ny);
Dauto1x(:,:,i)=reshape(sccsol.uncertainty(:,23),Nx,Ny);
Dauto1y(:,:,i)=reshape(sccsol.uncertainty(:,24),Nx,Ny);
Dauto2x(:,:,i)=reshape(sccsol.uncertainty(:,27),Nx,Ny);
Dauto2y(:,:,i)=reshape(sccsol.uncertainty(:,28),Nx,Ny);
D1x(:,:,i)=reshape(sccsol.uncertainty(:,31),Nx,Ny);
D1y(:,:,i)=reshape(sccsol.uncertainty(:,32),Nx,Ny);
end

% keyboard;
end

%% Save velocity, error and uncertainties
% save([savematdir,'2003BSOM.mat'],'err_sccu','err_sccv','IXX1','IYY1','IXX2','IYY2');

%% Postprocess
%%
%convert to vectors
err_sccu=err_sccu(:);
err_sccv=err_sccv(:);
noabserrx=noabserrx(:);
noabserry=noabserry(:);

MI=MI(:);
Autod=Autod(:);
IXX1=Ixx(:)./sqrt(MI);
IYY1=Iyy(:)./sqrt(MI);
IXX2=IXXd(:);
IYY2=IYYd(:);

Gpx=Gpx(:);
Gpy=Gpy(:);
Spx=Spx(:);
Spy=Spy(:);

mewx=abs(Spx-Gpx);
mewy=abs(Spy-Gpy);

IXX1=sqrt(IXX1.^2+mewx.^2);
IYY1=sqrt(IYY1.^2+mewy.^2);


if dstore==1
   DCCx=DCCx(:);
   DCCy=DCCy(:);
   DGccx=DGccx(:);
   DGccy=DGccy(:);
   Dauto1x=Dauto1x(:);
   Dauto1y=Dauto1y(:);
   Dauto2x=Dauto2x(:);
   Dauto2y=Dauto2y(:);
   D1x=D1x(:);
   D1y=D1y(:);
end
%%
%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)=[];

% upprx(inx,:)=[];uprmsx(inx,:)=[];upcex(inx,:)=[];uenx(inx,:)=[];umix(inx,:)=[];
% uppry(iny,:)=[];uprmsy(iny,:)=[];upcey(iny,:)=[];ueny(iny,:)=[];umiy(iny,:)=[];
%%
%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)]

noabserrx(inx2)=[];
noabserry(iny2)=[];
noabserrx(inx3)=[];
noabserry(iny3)=[];

%%
thresh=1;
%invalid vector removal
inx1=(find(err_sccu(:)>thresh));
iny1=(find(err_sccv(:)>thresh));
%keyboard;
IXX1(inx1)=[];
IYY1(iny1)=[];
IXX2(inx1)=[];
IYY2(iny1)=[];
err_sccu(inx1)=[];
err_sccv(iny1)=[];
[length(inx1) length(iny1)]
% upprx(inx1,:)=[];uprmsx(inx1,:)=[];upcex(inx1,:)=[];uenx(inx1,:)=[];umix(inx1,:)=[];
% uppry(iny1,:)=[];uprmsy(iny1,:)=[];upcey(iny1,:)=[];ueny(iny1,:)=[];umiy(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)=[];
%%
ll1=-0.4;lu1=0.4;
ll2=-0.2;lu2=0.2;
inx5=find(noabserrx>lu1 | noabserrx<ll1);
iny5=find(noabserry>lu2 | noabserry<ll2);
noabserrx(inx5)=[];
noabserry(iny5)=[];

[mu1,sig1]=normfit(noabserrx);
[mu2,sig2]=normfit(noabserry);
[mu1 sig1 mu2 sig2]

noabserrx1=noabserrx-mu1;
noabserry1=noabserry-mu2;

% % IXX1(inx5)=[];
% % IYY1(iny5)=[];
% % IXX2(inx5)=[];
% % IYY2(iny5)=[];
% % err_sccu(inx5)=[];
% % err_sccv(iny5)=[];

% err_sccu=abs(noabserrx1);
% err_sccv=abs(noabserry1);

figure;hist(noabserrx,60);
figure;hist(noabserry,60);
[mean(noabserrx1) mean(noabserry1)]
% keyboard;
%% Coverage Calculation
N=r*s*t;

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

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

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

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

[covx covx2]
[covy covy2]
%% Plots
% %Plot for deform
lw=3;
fs=18;
vec=linspace(0,0.1,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:8e4;%max(Nrdx);
% ll2=1:100:8e4;%max(Nrdy);

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


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

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

mixx2=rms(IXX2(IXX2<=0.1)).*ones(length(ll),1);
miyy2=rms(IYY2(IYY2<=0.1)).*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;
ylabel('No. of Data Points');
xlabel('pixels');
title(sprintf(['Ixx=',num2str(covx),'  IMx=',num2str(covx2)]));
if printfig==1
    print(gcf,'-dpng',[outdir,'rmsUhist.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;
ylabel('No. of Data Points');xlabel('pixels');
title(sprintf(['Iyy=',num2str(covy),'  IMy=',num2str(covy2)]));
if printfig==1
    print(gcf,'-dpng',[outdir,'rmsVhist.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');



Nbin=10;

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

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('Covergae 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']);