Skip to content
Permalink
a296dcf2a6
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

Moment_of_Correlation_2DPIV_uncertainty/postprocess_codes/Jet_data_postprocess.m

Go to file
 
 
Cannot retrieve contributors at this time
839 lines (710 sloc) 29.1 KB
Raw Blame
% Comapre Prana and Davis Solutions for experimental Jet data
clear all
strcase='cropped';
printfig=0;
plotfig=1;
savematdir='Z:\Planar_Uncertainty_work\Results\Matfiles\06_18_2015\';
N=496;
truesoldir='Z:\Jet_Uncertainty_Data\HDR_data_F001\';
% davisdir='Z:\Planar_Uncertainty_work\uncertainty_tests_Davis\uncertainty_test_jetflow\MS_images_F001_cropped_x_303_y_188_frame_03_501\TR_PIV_MP(3x16x16_75%ov)_01\'; % without vector pre and post processing
% davisdir='Z:\Planar_Uncertainty_work\uncertainty_tests_Davis\uncertainty_test_jetflow\MS_images_F001_cropped_x_303_y_188_frame_03_501\TR_PIV_MP(1x32x32_87%ov)(3x16x16_75%ov)_nopostprocess\'; % without vector pre and post processing
davisdir='Z:\Planar_Uncertainty_work\uncertainty_tests_Davis\uncertainty_test_jetflow\MS_images_F001_cropped_x_303to402_y_188to347_frame_03_501\TR_PIV_MP(1x32x32_75%ov)(3x16x16_75%ov)_postprocess\'; % without vector pre and post processing
% soldir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\';
% soldir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\newixx\';
% soldir2='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\Jetdata\';
% soldir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Jetdata\';
% soldir2='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\Jetdata\ws32IM\';
soldir='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\Jetdata\';
soldir2='Z:\Planar_Uncertainty_work\Results\10_22_2015_different_processing_windows\Jetdata\';
% soldir2='Z:\Planar_Uncertainty_work\Results\Experimental_Jet\unsteady_inviscid_core\New_processing_test_compare_prana_and_Davis\cropped_x_303_y_188\set3_500_firstpass64\';% Image Matching
soldircs='Z:\Planar_Uncertainty_work\Results\Experimental_Jet\unsteady_inviscid_core\New_processing_test_compare_prana_and_Davis\cropped_x_303_y_188\Correlation_Statistics\with_filter\'; % Corr stat
% outputdir='Z:\Planar_Uncertainty_work\Results\09_23_2015_final_processing\plots\Jetdata\newixx2\';
outputdir='Z:\Planar_Uncertainty_work\Results\10_14_2015_deform_whittaker_blackman\plots\Jetdata\ws32IM\';
biasfact=0;
if ~exist(outputdir,'dir')
mkdir(outputdir);
end
addpath Y:\Projects\TurbulentFlame\analysis\src\readimx_v2.0_win64\readimx;
basename1='PIV_scc_jet_cropped_pass4_';
basename2='B';
basename3='PIV_scc_jet_cropped_pass4_';
% c1=4;c2=26;%303+14-1 303+102-1
% r1=4; r2=36; %188+14-1 188+142-1
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;
Sx=40;Sy=25;
for i=1:N
truesol=importdata([truesoldir,'B',num2str(i+2,'%05.0f'),'.dat']);
sccsol=load([soldir,basename1,num2str(i+2,'%05.0f'),'.mat']);
sccsolcs=load([soldircs,basename3,num2str(i+2,'%05.0f'),'.mat']);
sccsol2=load([soldir2,basename3,num2str(i+2,'%05.0f'),'.mat']);
%% Davis Result
M=readimx([davisdir,basename2,num2str(i,'%05.0f'),'.vc7']);
if i==1
[D] = create2DVec(M.Frames{1,1});
Xd=D.X;
Yd=D.Y;
Xd=Xd';
Yd=Yd';
if strcmpi(strcase,'whole')
Xd=Xd(r1:r2,c1:c2);
Yd=Yd(r1:r2,c1:c2);
elseif strcmpi(strcase,'cropped')
Xd=Xd(r1:r2,c1:c2)+302;
Yd=Yd(r1:r2,c1:c2)+187;
end
clear D;
end
Udtemp=double(cell2mat(M.Frames{1,1}.Components{1,1}.Planes)');
Vdtemp=double(cell2mat(M.Frames{1,1}.Components{2,1}.Planes)');
Ud(:,:,i)=Udtemp(r1:r2,c1:c2);
Vd(:,:,i)=Vdtemp(r1:r2,c1:c2);
%% Get Correlation Statistics estimates
Unrxtemp=cell2mat(M.Frames{1,1}.Components{22,1}.Planes)';
Unrytemp=cell2mat(M.Frames{1,1}.Components{23,1}.Planes)';
Unbxtemp=cell2mat(M.Frames{1,1}.Components{25,1}.Planes)';
Unbytemp=cell2mat(M.Frames{1,1}.Components{26,1}.Planes)';
Untrtemp=cell2mat(M.Frames{1,1}.Components{21,1}.Planes)';
Untbtemp=cell2mat(M.Frames{1,1}.Components{24,1}.Planes)';
Unrx(:,:,i)=Unrxtemp(r1:r2,c1:c2);
Unry(:,:,i)=Unrytemp(r1:r2,c1:c2);
Unbx(:,:,i)=Unbxtemp(r1:r2,c1:c2);
Unby(:,:,i)=Unbytemp(r1:r2,c1:c2);
Untr(:,:,i)=Untrtemp(r1:r2,c1:c2);
Untb(:,:,i)=Untbtemp(r1:r2,c1:c2);
%% 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=6:2:71;
% tsy=5:2:50;
tsx=8:2:69;
tsy=7:2:43;
X=X(tsx,tsy);
Y=Y(tsx,tsy);
Xh=(Xd-105)./(10.2*7.2);
Yh=-((Yd-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);
%% 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);
%Scc peak location
Spx=reshape(sccsol.uncertainty(:,17),Sx,Sy);
Spx=Spx(end:-1:1,:);
Spx=biasfact*Spx(r1:r2,c1:c2);
Spy=reshape(sccsol.uncertainty(:,18),Sx,Sy);
Spy=Spy(end:-1:1,:);
Spy=biasfact*Spy(r1:r2,c1:c2);
%10:89,22:101;
%Bias error
mewx=abs(Spx-Gpx);
mewy=abs(Spy-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);
% IXX1(:,:,i)=Ixx;
% IYY1(:,:,i)=Iyy;
%% 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);
%% Get prana correlation statistics estimates
% tcsx=reshape(sccsolcs.uncertainty(:,25),Sx,Sy);
% tcsx=tcsx(end:-1:1,:);
% Uncsx(:,:,i)=tcsx(r1:r2,c1:c2);
%
% tcsy=reshape(sccsolcs.uncertainty(:,26),Sx,Sy);
% tcsy=tcsy(end:-1:1,:);
% Uncsy(:,:,i)=tcsy(r1:r2,c1:c2);
%% Error
error_Us(:,:,i)=(Utrue(:,:,i)-Us(:,:,i));
error_Vs(:,:,i)=(Vtrue(:,:,i)-Vs(:,:,i));
error_Ud(:,:,i)=(Utrue(:,:,i)-Ud(:,:,i));
error_Vd(:,:,i)=(Vtrue(:,:,i)-Vd(:,:,i));
UdUs(:,:,i)=(Ud(:,:,i)-Us(:,:,i));
VdVs(:,:,i)=(Vd(:,:,i)-Vs(:,:,i));
%% 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,Ud(:,:,i));axis image;shading interp;caxis([0 4]);colorbar;title('Udavis');
% % print(gcf,'-dpng',[outputdir,'Udavis_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;pcolor(Xh,Yh,Vd(:,:,i));axis image;shading interp;caxis([-0.6 0.6]);colorbar;title('Vdavis');
% % print(gcf,'-dpng',[outputdir,'Vdavis_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');
% % quiver(Xh,Yh,Ud(:,:,i),Vd(:,:,i),1,'color','blue');
% % 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_V(:,:,i));axis image;shading interp;caxis([-0.2 0.2]);title('Vtrue');colorbar;
% subplot(1,2,1);pcolor(Xh,Yh,error_U(:,:,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;
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');
% keyboard;
%%
teUs=error_Us(:);
% teUs(abs(teUs)>0.5)=[];
teUd=error_Ud(:);
% teUd(abs(teUd)>0.5)=[];
teUsd=UdUs(:);
% teUsd(abs(teUsd)>0.5)=[];
teVs=error_Vs(:);
% teVs(abs(teVs)>0.5)=[];
teVd=error_Vd(:);
% teVd(abs(teVd)>0.5)=[];
teVsd=VdVs(:);
% teVsd(abs(teVsd)>0.5)=[];
vec=linspace(-0.5,0.5,50);
Neus=histc(teUs,vec);
Neud=histc(teUd,vec);
Neusd=histc(teUsd,vec);
Nevs=histc(teVs,vec);
Nevd=histc(teVd,vec);
Nevsd=histc(teVsd,vec);
figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
hold on;
plot(vec,Neus,'b-');
plot(vec,Neud,'r-');
% plot(vec,Neusd,'k-');
% legend('Utrue-Uprana','Utrue-Udavis','Udavis-Uprana');
legend('Utrue-Uprana','Utrue-Udavis');
title('U error histogram');
hold off;
grid on;
% keyboard;
if printfig==1
print(gcf,'-dpng',[outputdir,'UerrorHist2','.png']);
end
figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
hold on;
plot(vec,Nevs,'b-');
plot(vec,Nevd,'r-');
% plot(vec,Nevsd,'k-');
% legend('Vtrue-Vprana','Vtrue-Vdavis','Vdavis-Vprana');
legend('Vtrue-Vprana','Vtrue-Vdavis');
title('V error histogram');
hold off;
grid on;
if printfig==1
print(gcf,'-dpng',[outputdir,'VerrorHist2','.png']);
end
% keyboard;
%% Instantaneous U Velocity profile
xloc=16;yloc=14; %(378,254, 376,253)
% xloc=21;yloc=14; %(378,254, 376,253)
tstamp=10;
Uinsttrue=squeeze(Utrue(:,xloc,tstamp));
Uinstsol=squeeze(Us(:,xloc,tstamp));
Uinstdav=squeeze(Ud(:,xloc,tstamp));
figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
plot(Yh(:,1),Uinsttrue,'k-');hold on;
plot(Yh(:,1),Uinstsol,'b-');
plot(Yh(:,1),Uinstdav,'r-');hold off;
axis([-1.4 1.4 -0.4 4.2]);
xlabel('Yh');
ylabel('U');
title(['Instantaneous U Velocity at t=',num2str(0.1*tstamp,'%02.1f'),' ms']);
legend('True','prana','davis');
box on;
% set(gcf,'units','normalized','outerposition',[0 0 1 1]);
if printfig==1
print(gcf,'-dpng',[outputdir,'InstantU.png']);
end
%% velocity time series
UHDts=squeeze(Utrue(yloc,xloc,:));
UMSts=squeeze(Us(yloc,xloc,:));
UDVts=squeeze(Ud(yloc,xloc,:));
ts=0:0.1:0.1*length(UHDts)-0.1;
ts=ts';
figure;set(gca,'FontSize',fs);
hold on;
plot(ts,UHDts,'k-','LineWidth',3);
plot(ts,UMSts,'b-','LineWidth',1.5);
plot(ts,UDVts,'r-','LineWidth',1.5);
axis([0 40 2.95 4.05]);
hold off;
xlabel('t (ms)');
ylabel('U (px)');
title(['Velocity time series at X=',num2str(X(1,xloc),'%02.1f'),'Y=',num2str(Y(yloc,1),'%02.1f')]);
legend('True','Prana','Davis');
set(gcf,'units','normalized','outerposition',[0 0 1 1]);
box on;
if printfig==1
print(gcf,'-dpng',[outputdir,'Velocity time series.png']);
end
%% Error Time Series
eUts=squeeze(error_Us(yloc,xloc,:));
eUdts=squeeze(error_Ud(yloc,xloc,:));
eVts=squeeze(error_Vs(yloc,xloc,:));
eMts=sqrt(eUts.^2+eVts.^2);
meUts=mean(eUts).*ones(length(eUts));
meUdts=mean(eUdts).*ones(length(eUdts));
figure;set(gca,'FontSize',fs);
plot(ts,eUts,'b-');hold on;
plot(ts,meUts,'b--','Linewidth',2);hold off;
axis([0 max(ts(:)) -0.19 0.19]);
xlabel('t (ms)');
ylabel('e_U prana(px)');
box on;
st1='Error U time series (prana)';
st4='Error U time series (Davis)';
st2=['At X=',num2str(X(1,xloc),'%02.1f'),',Y=',num2str(Y(yloc,1),'%02.1f')];
st3=['mean error=',num2str(meUts(1),'%02.3f')];
st5=['mean error=',num2str(meUdts(1),'%02.3f')];
title([st1,'\newline',st2,'\newline',st3]);
if printfig==1
print(gcf,'-dpng',[outputdir,'Error U (Prana) time series.png']);
end
figure;set(gca,'FontSize',fs);
plot(ts,eUdts,'r-');hold on;
plot(ts,meUdts,'r--','Linewidth',2);hold off;
axis([0 max(ts(:)) -0.19 0.19]);box on;
xlabel('t (ms)');
ylabel('e_U davis(px)');
title([st4,'\newline',st2,'\newline',st5]);
if printfig==1
print(gcf,'-dpng',[outputdir,'Error U (Davis) time series.png']);
end
%% Instantaneous Velocity and Error profile
etemp1=squeeze(error_Us(:,xloc,:));
rmserrU=rms(etemp1,2);
etemp5=squeeze(error_Ud(:,xloc,:));
rmserrUd=rms(etemp5,2);
etemp2=squeeze(IXX1(:,xloc,:));
rmserrIxx=rms(etemp2,2);
etemp3=squeeze(Imx(:,xloc,:));
rmserrImx=rms(etemp3,2);
etemp6=squeeze(Unrx(:,xloc,:));
rmserrIxxd1=rms(etemp6,2);
etemp7=squeeze(Unbx(:,xloc,:));
rmserrIxxd2=rms(etemp7,2);
rmserrIxxd=sqrt(rmserrIxxd1.^2+rmserrIxxd2.^2);
% etemp8=squeeze(Uncsx(:,xloc,:));
% rmserrUncsx=rms(etemp8,2);
figure;set(gcf,'DefaultLineLineWidth',lw);set(gca,'FontSize',fs);
hold on;
plot(Yh(:,1),rmserrU,'r-');
plot(Yh(:,1),rmserrUd,'r--');
plot(Yh(:,1),rmserrIxx,'k-');
plot(Yh(:,1),rmserrImx,'c-');
plot(Yh(:,1),rmserrIxxd,'m--');
% plot(Yh(:,1),rmserrUncsx,'m-');
hold off;box on;
axis([-1.4 1.4 0 0.2]);
xlabel('Yh');
ylabel('RMS error, uncertainty');
title(['RMS error over time at X=',num2str(X(1,xloc),'%02.1f')]);
% legend('rms e_u (Prana)','rms e_u (Davis)','rms IXX (Prana)','rms PDx(Prana)','rms UCSx(Davis)','rms UCSx (Prana)');
set(gcf,'units','normalized','outerposition',[0 0 1 1]);
if printfig==1
print(gcf,'-dpng',[outputdir,'RMS error2.png']);
end
keyboard;
%%
if plotfig==1
error_U=abs(error_Us);
error_V=abs(error_Vs);
%% 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.5));
% iny1=(find(error_V(:)>0.5));
%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));
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.5;
vec=linspace(0,thresh,60);
% 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(error_U,vec);
Nrey=histc(error_V,vec);
Nrdx=histc(IXX1,vec);
Nrdy=histc(IYY1,vec);
Nrdx2=histc(Imx,vec);
Nrdy2=histc(Imy,vec);
% Npprxl=histc(upprx(:,2),vec2);
% Npprxh=histc(upprx(:,1),vec3);
%
% Nppryl=histc(uppry(:,2),vec2);
% Nppryh=histc(uppry(:,1),vec3);
%
% Nmixl=histc(umix(:,2),vec4);
% Nmixh=histc(umix(:,1),vec5);
%
% Nmiyl=histc(umiy(:,2),vec4);
% Nmiyh=histc(umiy(:,1),vec5);
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=median(IXX1(IXX1<=thresh)).*ones(length(ll),1);
% miyy1=median(IYY1(IYY1<=thresh)).*ones(length(ll2),1);
%
% mixx2=median(Imx(Imx<=thresh)).*ones(length(ll),1);
% miyy2=median(Imy(Imy<=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--');
% plot(vec2,Npprxl,'color',rgb('DarkRed'));
% plot(vec3,Npprxh,'color',rgb('DarkRed'));
% plot(vec4,Nmixl,'color',rgb('orange'));
% plot(vec5,Nmixh,'color',rgb('orange'));
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)');
%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--');
% 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)');
%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)]
%%
plotflag=1;
if plotflag==1;
IXX2=Imx;
IYY2=Imy;
err_sccu=error_U;
err_sccv=error_V;
Nbin=15;
maxixx=thresh;
minixx=0;
vu=linspace(minixx,maxixx,Nbin);
maxiyy=thresh;
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',[outputdir,'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',[outputdir,'rms error in binned uncertainty.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',[outputdir,'covg in each bin.png']);
end
end
end