stagnation_flow_postprocess.m

% Testing IXX for STagnation Flow
clear all;
%Load true solution:-
Ut=load('Z:\Planar_Uncertainty_work\codes\matfiles\Utrue.mat');
Vt=load('Z:\Planar_Uncertainty_work\codes\matfiles\Vtrue.mat');
Utrue=Ut.Ue;
Vtrue=Vt.Ve;



Un_fit=load('Z:\Planar_Uncertainty_work\codes\matfiles\fit_un.mat');
unfitx=Un_fit.u_pUPe1;
unfity=Un_fit.u_pVPe1;
% figure;imagesc(Ut1);colorbar;caxis([-0.06 0.06]);
% figure;imagesc(Vt1);colorbar;caxis([0 3.5]);

istring1=sprintf(['%%s%%0%0.0fd.','mat'],6);

% dirsol='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\stagnation_flow\PIV_scc_finalstag_pass4_';
% dirsol2='Z:\Planar_Uncertainty_work\Results\07_20_2015_scaling_tests\NewIM\stagnation_flow\PIV_scc_finalstag_pass5_';

% dirsol='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\stagnation_flow\PIV_scc_finalstag_pass4_';
% dirsol='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\stagnation_flow\newixx\PIV_scc_finalstag_pass4_';
% dirsol2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\stagnation_flow\PIV_scc_finalstag_pass5_';
dirsol='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\stagnation_flow\PIV_scc_finalstag_pass4_';
dirsol2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\stagnation_flow\ws32IM\PIV_scc_finalstag_pass4_';


% outdir1='Z:\Planar_Uncertainty_work\Results\Ixx_with_cclsq_dia\plots\Stagnation_flow\';
% outdir1='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\plots\stagnation_flow\';
% outdir1='Z:\Planar_Uncertainty_work\Results\07_07_2015_tests_new_scaling\With_zero_center_gauss\plots\stagnation_flow\';
% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\stagnation_flow\';
% outdir1='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\stagnation_flow\newixx\';
outdir1='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\stagnation_flow\ws32IM\';

biasfact=0;

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

t=99;

Sx=64; 
Sy=80;

for i=1:t
    sccsol=load(sprintf(istring1,dirsol,i));
    sccsol2=load(sprintf(istring1,dirsol2,i));
    
    %Solution
    Us(:,:,i)=sccsol.U;
    Vs(:,:,i)=sccsol.V;
    
    %Error 
    err_sccu(:,:,i)=abs(Utrue-Us(:,:,i));
    err_sccv(:,:,i)=abs(Vtrue-Vs(:,:,i));
    
    % Image Matching
    Imx=reshape(sccsol2.uncertainty(:,9),Sx,Sy);
    Imy=reshape(sccsol2.uncertainty(:,10),Sx,Sy);
    
    MI(:,:,i)=reshape(sccsol.uncertainty(:,5),Sx,Sy);
    autod=reshape(sccsol.uncertainty(:,6),Sx,Sy);
    Ixx=reshape(sccsol.uncertainty(:,7),Sx,Sy);
    Iyy=reshape(sccsol.uncertainty(:,8),Sx,Sy);
    
   % correct gradient (if desired)
    [Udiff]=socdiff(Us(:,:,i),16,1);
    [Vdiff]=socdiff(Vs(:,:,i),16,2);
    Ixx = real(sqrt(Ixx.^2 - (autod.^2/16).*(Udiff).^2 ));
    Iyy = real(sqrt(Iyy.^2 - (autod.^2/16).*(Vdiff).^2 ));
    
   %Gcc peak location  
    Gpx=reshape(sccsol.uncertainty(:,15),Sx,Sy);
    Gpy=reshape(sccsol.uncertainty(:,16),Sx,Sy);
   
    %Scc peak location
    Spx=biasfact*reshape(sccsol.uncertainty(:,17),Sx,Sy);
    Spy=biasfact*reshape(sccsol.uncertainty(:,18),Sx,Sy);
    
    %Bias error
    mewx=abs(Spx-Gpx);
    mewy=abs(Spy-Gpy);

    
    %Scaling for ixx
    Ixx=Ixx./sqrt(MI(:,:,i));
    Iyy=Iyy./sqrt(MI(:,:,i));
        
    % Final second order moment uncertainty
    IXX1(:,:,i)=sqrt(Ixx.^2+mewx.^2+unfitx.^2);
    IYY1(:,:,i)=sqrt(Iyy.^2+mewy.^2+unfity.^2);
    % Image Matching
    IXX2(:,:,i)=sqrt(Imx.^2+unfitx.^2);
    IYY2(:,:,i)=sqrt(Imy.^2+unfity.^2);
    
end
%% convert to vectors
IXX1=IXX1(:);
IYY1=IYY1(:);
IXX2=IXX2(:);
IYY2=IYY2(:);
err_sccu=err_sccu(:);
err_sccv=err_sccv(:);
MI=MI(:);
%%
%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)=[];
%%
%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
inx1=(find(err_sccu(:)>1));
iny1=(find(err_sccv(:)>1));
%keyboard;
IXX1(inx1)=[];
IYY1(iny1)=[];
IXX2(inx1)=[];
IYY2(iny1)=[];
err_sccu(inx1)=[];
err_sccv(iny1)=[];

%%

gmx=(find(err_sccu(:)<=1));
gmy=(find(err_sccv(:)<=1));
ngmx=length(gmx);
ngmy=length(gmy);
% 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]% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
lw=3;
fs=18;
 %% Plot error and uncertainty histogram
thresh=0.3;
 
vec=linspace(0,thresh,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:max(Nrex);
ll2=1:100:max(Nrey);


sigu=rms(err_sccu(err_sccu<=thresh)).*ones(length(ll),1);
sigv=rms(err_sccv(err_sccv<=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(IXX2(IXX2<=thresh)).*ones(length(ll),1);
miyy2=rms(IYY2(IYY2<=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 thresh 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 thresh 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;
    
    Nbin=10;
    
    maxixx=0.15;
    minixx=0.02;
    vu=linspace(minixx,maxixx,Nbin);
    maxiyy=0.3;
    miniyy=0.04;
    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,'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