Jet_data_postprocess_pranaonly.m

% Comapre Prana and Davis Solutions for experimental Jet data
clear all
strcase='cropped';
printfig=0;
N=496;
% Previous
%X 316 to 404 (88); Y 201 to 329 (128)
% Now
% X 320 to 392 (72); Y 205 to 325 (120)
c1=5;c2=23;%303+14-1  303+102-1
r1=5; r2=35; %188+14-1  188+142-1
lw=3;fs=20;
I=78;J=70;
% Sx=66;Sy=46;
Sx=40;Sy=30;

% INPUT DIRECTORY
% soldir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\PIV_scc_jet_cropped_pass4_';
% soldir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\newixx\PIV_scc_jet_cropped_pass4_';
% soldir2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\PIV_scc_jet_cropped_pass5_';

% soldir='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\Jetdata\PIV_scc_jet_cropped_pass4_';
% soldir2='Z:\Planar_Uncertainty_work\Results\10_05_2015_Newdeform_NewIM\Jetdata\PIV_scc_jet_cropped_pass5_';
soldir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Jetdata\PIV_scc_jet_cropped_pass4_';
soldir2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Jetdata\ws32IM\PIV_scc_jet_cropped_pass4_';


% Output Directory
savematdir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Matfiles\Jetdata\pranaonly\';
% outputdir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\Jetdata\pranaonly\newixx\';
outputdir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\Jetdata\pranaonly\ws32IM\';
biasfact=0;

if ~exist(savematdir,'dir')
    mkdir(savematdir);
end
if ~exist(outputdir,'dir')
    mkdir(outputdir);
end

% True solution directory
truesoldir='Z:\Jet_Uncertainty_Data\HDR_data_F001\B';


for i=1:N
    
    truesol=importdata([truesoldir,num2str(i+2,'%05.0f'),'.dat']);
    sccsol=load([soldir,num2str(i+2,'%05.0f'),'.mat']);
    sccsol2=load([soldir2,num2str(i+2,'%05.0f'),'.mat']);
%    keyboard;
    %% True solution
    if i==1
        X=truesol.data(:,1);
        Y=truesol.data(:,2);
        flag=truesol.data(:,6);
        X=reshape(X,J,I);
        Y=reshape(Y,J,I);
        X=X';Y=Y';
        tsx=8:2:69;
        tsy=7:2:43;
        X=X(tsx,tsy);
        Y=Y(tsx,tsy);
        
        Xh=(X-105)./(10.2*7.2);
        Yh=-((Y-254)./(10.2*7.2));
    end
    
    U=truesol.data(:,3);
    V=truesol.data(:,4);
    Ut=reshape(U,J,I)';
    Vt=reshape(V,J,I)';
    
    
    Utemp=Ut(tsx,tsy);
    Vtemp=Vt(tsx,tsy);
    
    
    
    
    Utrue(:,:,i)=Utemp;
    Vtrue(:,:,i)=Vtemp;
    
    %% prana Solution
    
    if i==1
        Xs=sccsol.X+302;
        Ys=sccsol.Y+187;
        Xs=Xs(r1:r2,c1:c2);
        Ys=Ys(r1:r2,c1:c2);
    end
    Ustemp=sccsol.U(end:-1:1,:);
    Vstemp=-sccsol.V(end:-1:1,:);
    Us(:,:,i)=Ustemp(r1:r2,c1:c2);
    Vs(:,:,i)=Vstemp(r1:r2,c1:c2);
    
%     keyboard;
    %% Get IXX/IYY estimates
    
    % Storing Ixx and other parameters.
    tppr=reshape(sccsol.uncertainty(:,1),Sx,Sy);
    tppr=tppr(end:-1:1,:);
    PPR(:,:,i)=tppr(r1:r2,c1:c2);
    
    tmi=reshape(sccsol.uncertainty(:,5),Sx,Sy);
    tmi=tmi(end:-1:1,:);
    MI(:,:,i)=tmi(r1:r2,c1:c2);
    
    Autod=reshape(sccsol.uncertainty(:,6),Sx,Sy);
    Autod=Autod(end:-1:1,:);
    Autod=Autod(r1:r2,c1:c2);
    
    Ixx=reshape(sccsol.uncertainty(:,7),Sx,Sy);
    Ixx=Ixx(end:-1:1,:);
    Ixx=Ixx(r1:r2,c1:c2);
    
    Iyy=reshape(sccsol.uncertainty(:,8),Sx,Sy);
    Iyy=Iyy(end:-1:1,:);
    Iyy=Iyy(r1:r2,c1:c2);
    
    %Gcc peak location  
    Gpx=reshape(sccsol.uncertainty(:,15),Sx,Sy);
    Gpx=Gpx(end:-1:1,:);
    Gpx=Gpx(r1:r2,c1:c2);
    
    Gpy=reshape(sccsol.uncertainty(:,16),Sx,Sy);
    Gpy=Gpy(end:-1:1,:);
    Gpy=Gpy(r1:r2,c1:c2);
            
    %Bias error
    mewx=abs(Gpx);
    mewy=abs(Gpy);
    % correct for gradient
    [Udiff]=socdiff(Us(:,:,i),2,1);
    [Vdiff]=socdiff(Vs(:,:,i),2,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);

   
    %% Get IMx/IMy estimates
    % Image Matching uncertainty
    timx=reshape(sccsol2.uncertainty(:,9),Sx,Sy);
    timx=timx(end:-1:1,:);
    Imx(:,:,i)=timx(r1:r2,c1:c2);
    
    timy=reshape(sccsol2.uncertainty(:,10),Sx,Sy);
    timy=timy(end:-1:1,:);
    Imy(:,:,i)=timy(r1:r2,c1:c2);
    
    tnim=reshape(sccsol2.uncertainty(:,11),Sx,Sy);
    tnim=tnim(end:-1:1,:);
    Nim(:,:,i)=tnim(r1:r2,c1:c2);
    
    
    %% Error
    error_Us(:,:,i)=(Utrue(:,:,i)-Us(:,:,i));
    error_Vs(:,:,i)=(Vtrue(:,:,i)-Vs(:,:,i));
    
    %%
%     figure(1);
%     subplot(2,3,1);imagesc(abs(error_Us(:,:,i)));colorbar;caxis([0 0.2]);
%     subplot(2,3,2);imagesc(IXX1(:,:,i));colorbar;caxis([0 0.3]);
%     subplot(2,3,3);imagesc(Imx(:,:,i));colorbar;caxis([0 0.3]);
%     subplot(2,3,4);imagesc(abs(error_Vs(:,:,i)));colorbar;caxis([0 0.2]);
%     subplot(2,3,5);imagesc(IYY1(:,:,i));colorbar;caxis([0 0.3]);
%     subplot(2,3,6);imagesc(Imy(:,:,i));colorbar;caxis([0 0.3]);
%     pause(0.2);
    %% Plot Velcity distributions
    %{
% %     if printfig==1
%         if i==1%>=1 || i<=10
            figure;pcolor(Xh,Yh,Utrue(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Utrue');
%             print(gcf,'-dpng',[outputdir,'Utrue_frame_',num2str(i+2,'%05.0f'),'.png']);
            figure;pcolor(Xh,Yh,Us(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Uprana');
%             print(gcf,'-dpng',[outputdir,'Uprana_frame_',num2str(i+2,'%05.0f'),'.png']);
            
            figure;pcolor(Xh,Yh,Vtrue(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vtrue');
%             print(gcf,'-dpng',[outputdir,'Vtrue_frame_',num2str(i+2,'%05.0f'),'.png']);
            figure;pcolor(Xh,Yh,Vs(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vprana');
%             print(gcf,'-dpng',[outputdir,'Vprana_frame_',num2str(i+2,'%05.0f'),'.png']);
             
            figure;%set(gcf,'DefaultLineLineWidth',1);
            set(gca,'FontSize',10);hold on;
            quiver(Xh,Yh,Utrue(:,:,i),Vtrue(:,:,i),1,'color','red');axis image; 
            quiver(Xh,Yh,Us(:,:,i),Vs(:,:,i),1,'color','black');
            hold off;
            box on;
            xlabel('x/h');ylabel('y/h');
            title(['Velocity field at frame',num2str(i+2,'%05.0f')]);
%             legend('True','prana','davis');
%             print(gcf,'-dpng',[outputdir,'Vectorfield_frame2_',num2str(i+2,'%05.0f'),'.png'],'-r1200');
            
%             keyboard;

%             pause(0.1);
%             close all;
%             pause(0.2);
%         end
%     end
%         if i==1
        figure;
        subplot(1,2,2);pcolor(Xh,Yh,error_Vs(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Vtrue');colorbar;
        subplot(1,2,1);pcolor(Xh,Yh,error_Us(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Utrue');colorbar;
%         if printfig==1
%             set(gcf,'units','normalized','outerposition',[0 0 1 1]);
%             print(gcf,'-dpng',[outputdir,'Error_frame_',num2str(i+2,'%05.0f'),'.png'],'-r300');
%         end
%     end

    keyboard;
    %}
    keyboard;
    if rem(i,100)==0
        i
    end
end
%% Save error and uncertainties
% save([savematdir,'Jet_uncertainty_data_x_303_y_188_firstpass64_cclsq.mat'],'error_Us','error_Vs','error_Ud','error_Vd','IXX1','IYY1','Imx','Imy','Unrx','Unry','Unbx','Unby');

error_U=abs(error_Us);
error_V=abs(error_Vs);

Ixxtemp=IXX1;
Iyytemp=IYY1;

 figure(1);
 subplot(2,3,1);imagesc(mean(error_U,3));colorbar;caxis([0 0.2]);
 subplot(2,3,2);imagesc(mean(IXX1,3));colorbar;caxis([0 0.2]);
 subplot(2,3,3);imagesc(nanmean(Imx,3));colorbar;caxis([0 0.1]);
 subplot(2,3,4);imagesc(mean(error_V,3));colorbar;caxis([0 0.1]);
 subplot(2,3,5);imagesc(mean(IYY1,3));colorbar;caxis([0 0.1]);
 subplot(2,3,6);imagesc(nanmean(Imy,3));colorbar;caxis([0 0.1]);
%     pause(0.2);

if printfig==1
    print(gcf,'-dpng',[outputdir,'Mean_profiles.png']);
end

% keyboard;




%% convert to vectors
IXX1=IXX1(:);
IYY1=IYY1(:);
Imx=Imx(:);
Imy=Imy(:);
error_U=error_U(:);
error_V=error_V(:);
% IXX2=Imx;
% IYY2=Imy;
% 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)=[];
%%
%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
inx1=(find(error_U(:)>1));
iny1=(find(error_V(:)>1));
% inx1=(find(error_U(:)>0.4));
% iny1=(find(error_V(:)>0.4));
%keyboard;
IXX1(inx1)=[];
IYY1(iny1)=[];
Imx(inx1)=[];
Imy(iny1)=[];
error_U(inx1)=[];
error_V(iny1)=[];

%% Find Coverage

% gmx=(find(error_U(:)<=1));
% gmy=(find(error_V(:)<=1));
nx=length(error_U);
ny=length(error_V);

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

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

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

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

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

[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.4;

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);%max(Nrdx);
ll2=1:100:max(Nrey);%max(Nrdy);


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(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',[outputdir,'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',[outputdir,'Vhist2.png']);
end
[sigu(1) mixx1(1) mixx2(1)]
[sigv(1) miyy1(1) miyy2(1)]


IXX2=Imx;
IYY2=Imy;
err_sccu=error_U;
err_sccv=error_V;

Nbin=10;

maxixx=0.2;
minixx=0;
vu=linspace(minixx,maxixx,Nbin);
maxiyy=0.2;
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)=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(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',[outputdir,'rms error in binned uncertainty.png']);
end

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(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',[outputdir,'covg in each bin.png']);
end
%%
% vecerr=linspace(-0.4,0.4,60);
% Nex=histc(error_Us(:),vecerr);
% Ney=histc(error_Vs(:),vecerr);
% 
% figure;plot(vecerr,Nex,'k-');
% figure;plot(vecerr,Ney,'k-');
% 
% errUt=error_Us(:);
% errVt=error_Vs(:);
% errUt1=errUt(abs(errUt)<0.4);
% errVt1=errVt(abs(errVt)<0.4);
% 
% [m1,s1]=normfit(errUt1)
% [m2,s2]=normfit(errVt1)

%%
% [~,indx]=sort(IXX1,'ascend');

% figure;plot(IXX1(indx),err_sccu(indx),'o');
% axis([0 0.2 0 0.3]);

% % figure;plot(err_sccu,'r.');
% % hold on;
% % plot(IXX1,'b.');hold off;