Index: /trunk/src/series/civ2vel_3C.m =================================================================== --- /trunk/src/series/civ2vel_3C.m (revision 877) +++ /trunk/src/series/civ2vel_3C.m (revision 878) @@ -169,5 +169,5 @@ %% MAIN LOOP ON FIELDS warning off -if NbField<2 && length(filecell{:,1})>2 +if NbField<2 disp_uvmat('ERROR','you need at least 2 images to compute the mean position for the stereo.',checkrun) return @@ -285,11 +285,10 @@ end ZI=mean(ZItemp,3); %mean between two the two time step + Vtest=ZItemp(:,:,2)-ZItemp(:,:,1); [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view a [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b - %trouver z - % trouver les coordonnées px sur chaque image. - %A= + for iview=1:2 %% reading input file(s) @@ -322,7 +321,9 @@ Ua=griddata(X1,Y1,U1,Xa,Ya); Va=griddata(X1,Y1,V1,Xa,Ya); - A=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,YI,YI,ZI); - - + + [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X + [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ + + %remove wrong vector temp=find(Data{2}.FF==0); X2=Data{2}.X(temp); @@ -332,7 +333,13 @@ Ub=griddata(X2,Y2,U2,Xb,Yb); Vb=griddata(X2,Y2,V2,Xb,Yb); - B=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,YI,YI,ZI); + + [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X + [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~ + + + % System to solve S=ones(size(XI,1),size(XI,2),3); D=ones(size(XI,1),size(XI,2),3,3); + S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; @@ -349,5 +356,5 @@ for indj=1:size(XI,1) for indi=1:size(XI,2) - dxyz=squeeze(S(indj,indi,:))\squeeze(D(indj,indi,:,:)); + dxyz=(squeeze(D(indj,indi,:,:))*1000)\(squeeze(S(indj,indi,:))*1000); % solving... U(indj,indi)=dxyz(1); V(indj,indi)=dxyz(2); @@ -360,4 +367,7 @@ Error(:,:,3)=B(:,:,1,1).*U+B(:,:,1,2).*V+B(:,:,1,3).*W-Ub; Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb; + + + %% generating the name of the merged field @@ -396,5 +406,8 @@ MergeData.V=V/Dt; MergeData.W=W/Dt; - MergeData.Error=sqrt(sum(Error.*Error,3)); + + mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; + mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; + MergeData.Error=(sqrt(mfx^2+mfy^2)/4).*sqrt(sum(Error.*Error,3)); errormsg=struct2nc(OutputFile,MergeData);%save result file if isempty(errormsg) @@ -406,8 +419,9 @@ -function A=get_coeff(Calib,X,Y,x,y,z) +function [A]=get_coeff(Calib,X,Y,x,y,z) % compute A~ coefficients R=(Calib.R)';%rotation matrix T_z=Calib.Tx_Ty_Tz(3); T=R(7)*x+R(8)*y+R(9)*z+T_z; + A(:,:,1,1)=(R(1)-R(7)*X)./T; A(:,:,1,2)=(R(2)-R(8)*X)./T; @@ -417,8 +431,24 @@ A(:,:,2,3)=(R(6)-R(9)*Y)./T; - - - -function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) +function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U + +X1d=Xd-Ud/2; +X2d=Xd+Ud/2; +Y1d=Yd-Vd/2; +Y2d=Yd+Vd/2; + +X1=(X1d-Calib.Cx_Cy(1))./Calib.fx_fy(1).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X1d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y1d-Calib.Cx_Cy(2)).^2 ).^(-1); +X2=(X2d-Calib.Cx_Cy(1))./Calib.fx_fy(1).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X2d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y2d-Calib.Cx_Cy(2)).^2 ).^(-1); +Y1=(Y1d-Calib.Cx_Cy(2))./Calib.fx_fy(2).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X1d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y1d-Calib.Cx_Cy(2)).^2 ).^(-1); +Y2=(Y2d-Calib.Cx_Cy(2))./Calib.fx_fy(2).*(1 + Calib.kc.*Calib.fx_fy(1).^(-2).*(X2d-Calib.Cx_Cy(1)).^2 + Calib.kc.*Calib.fx_fy(2).^(-2).*(Y2d-Calib.Cx_Cy(2)).^2 ).^(-1); + +U=X2-X1; +V=Y2-Y1; +X=X1+U/2; +Y=Y1+V/2; + + + +function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) % get H from stereo data z=0; error=0; @@ -429,5 +459,5 @@ R=(Calib_A.R)'; x_a=xmid- u/2; -y_a=ymid- u/2; +y_a=ymid- v/2; z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3); Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a; @@ -445,6 +475,10 @@ %% second image +%loading shift angle + Calib_B=XmlData{2}.GeometryCalib; R=(Calib_B.R)'; + + x_b=xmid+ u/2; y_b=ymid+ v/2; @@ -464,5 +498,11 @@ %% result Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya); -error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den; +error=((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den; +% ex=-error.*(Dyb-Dya); +% ey=-error.*(Dxb-Dxa); + +% z1=-u./(Dxb-Dxa); +% z2=-v./(Dyb-Dya); +z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den; xnew(1,:)=Dxa.*z+x_a; @@ -470,8 +510,8 @@ ynew(1,:)=Dya.*z+y_a; ynew(2,:)=Dyb.*z+y_b; - -Xphy=mean(xnew,1); %on moyenne les 2 valeurs -Yphy=mean(ynew,1); -z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den; - - +Xphy=mean(xnew,1); +Yphy=mean(ynew,1); + + + + Index: /trunk/src/series/civ_series.m =================================================================== --- /trunk/src/series/civ_series.m (revision 877) +++ /trunk/src/series/civ_series.m (revision 878) @@ -324,4 +324,5 @@ %%%%% MAIN LOOP %%%%%% maskname='';% initiate the mask name +tic; for ifield=1:NbField if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series (checkrun=1) @@ -852,4 +853,5 @@ end end +disp(['ellapsed time ' num2str(toc) ' s']) % 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/ @@ -1024,14 +1026,11 @@ XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); - image1_crop_1=image1_crop; image1_crop=interp2(image1_crop,XIant,YIant); - image1_crop_2=image1_crop; image1_crop(isnan(image1_crop))=0; - image1_crop_3=image1_crop; xi=(1:mesh:size(image2_crop,2)); yi=(1:mesh:size(image2_crop,1))'; - image2_crop=interp2(image2_crop,xi,yi); + image2_crop=interp2(image2_crop,xi,yi,'*spline'); image2_crop(isnan(image2_crop))=0; - par_civ.CorrSmooth=3; + %par_civ.CorrSmooth=3;%%%%%%%%%%%%%%%%%%% %% end Index: /trunk/src/series/stereo_civ.m =================================================================== --- /trunk/src/series/stereo_civ.m (revision 877) +++ /trunk/src/series/stereo_civ.m (revision 878) @@ -1,3 +1,3 @@ -%'stereo_civ': determination of topography by image correlation of two stereo views +%'civ_series': PIV function activated by the general GUI series % --- call the sub-functions: % civ: PIV function itself @@ -282,5 +282,6 @@ %% Civ1 - + + % if Civ1 computation is requested if isfield (Param.ActionInput,'Civ1') @@ -297,4 +298,7 @@ end end + + + [A,Rangx,Rangy]=phys_ima(A,XmlData,1); [Npy,Npx]=size(A{1}); @@ -516,15 +520,13 @@ mask=imread(par_civ2.Mask); end - ibx2=ceil(par_civ2.CorrBoxSize(1)/2); - iby2=ceil(par_civ2.CorrBoxSize(2)/2); -% par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess -% par_civ2.SearchBoxSize(2)=2*iby2+9; +% ibx2=ceil(par_civ2.CorrBoxSize(1)/2); +% iby2=ceil(par_civ2.CorrBoxSize(2)/2); par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors if par_civ2.CheckDeformation - par_civ2.DUDX=DUDX./nbval; - par_civ2.DUDY=DUDY./nbval; - par_civ2.DVDX=DVDX./nbval; - par_civ2.DVDY=DVDY./nbval; + par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1); + par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1); + par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1); + par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); end % calculate velocity data (y and v in indices, reverse to y component) @@ -617,8 +619,4 @@ Data.Civ2_FF(ind_good)=FFres; Data.CivStage=Data.CivStage+1; - - - - end @@ -659,7 +657,5 @@ [GridX,GridY]=meshgrid(minix:par_civ3.Dx:maxix,miniy:par_civ3.Dy:maxiy); par_civ3.Grid(:,1)=reshape(GridX,[],1); - par_civ3.Grid(:,2)=reshape(GridY,[],1); - - + par_civ3.Grid(:,2)=reshape(GridY,[],1); end Shiftx=zeros(size(par_civ3.Grid,1),1);% shift expected from civ2 data @@ -696,15 +692,13 @@ mask=imread(par_civ3.Mask); end - ibx2=ceil(par_civ3.CorrBoxSize(1)/2); - iby2=ceil(par_civ3.CorrBoxSize(2)/2); -% par_civ3.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess -% par_civ3.SearchBoxSize(2)=2*iby2+9; +% ibx2=ceil(par_civ3.CorrBoxSize(1)/2); +% iby2=ceil(par_civ3.CorrBoxSize(2)/2); par_civ3.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; par_civ3.Grid=[par_civ3.Grid(nbval>=1,1)-par_civ3.SearchBoxShift(:,1)/2 par_civ3.Grid(nbval>=1,2)-par_civ3.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors if par_civ3.CheckDeformation - par_civ3.DUDX=DUDX./nbval; - par_civ3.DUDY=DUDY./nbval; - par_civ3.DVDX=DVDX./nbval; - par_civ3.DVDY=DVDY./nbval; + par_civ3.DUDX=DUDX(nbval>=1)./nbval(nbval>=1); + par_civ3.DUDY=DUDY(nbval>=1)./nbval(nbval>=1); + par_civ3.DVDX=DVDX(nbval>=1)./nbval(nbval>=1); + par_civ3.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); end % calculate velocity data (y and v in indices, reverse to y component) @@ -724,6 +718,6 @@ nbvar=numel(Data.ListVarName); - Data.ListVarName=[Data.ListVarName {'Civ3_X','Civ3_Y','Civ3_U','Civ3_V','Civ3_F','Civ3_C','Xphys','Yphys','Zphys','Civ3_E'}];% cell array containing the names of the fields to record - Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}]; + Data.ListVarName=[Data.ListVarName {'Civ3_X','Civ3_Y','Civ3_U','Civ3_V','Civ3_F','Civ3_C','Xphys','Yphys','Zphys'}];% cell array containing the names of the fields to record + Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}]; Data.VarAttribute{nbvar+1}.Role='coord_x'; Data.VarAttribute{nbvar+2}.Role='coord_y'; @@ -776,7 +770,7 @@ Data.Patch3_SubDomainSize=Param.ActionInput.Patch3.SubDomainSize; nbvar=length(Data.ListVarName); - Data.ListVarName=[Data.ListVarName {'Civ3_U_smooth','Civ3_V_smooth','Civ3_SubRange','Civ3_NbCentres','Civ3_Coord_tps','Civ3_U_tps','Civ3_V_tps','Xmid','Ymid','Uphys','Vphys'}]; + Data.ListVarName=[Data.ListVarName {'Civ3_U_smooth','Civ3_V_smooth','Civ3_SubRange','Civ3_NbCentres','Civ3_Coord_tps','Civ3_U_tps','Civ3_V_tps','Xmid','Ymid','Uphys','Vphys','Error'}]; Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3',{'nb_coord','nb_bounds','nb_subdomain_3'},{'nb_subdomain_3'},... - {'nb_tps_3','nb_coord','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}]; + {'nb_tps_3','nb_coord','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},{'nb_tps_3','nb_subdomain_3'},'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}]; Data.VarAttribute{nbvar+1}.Role='vector_x'; @@ -805,6 +799,6 @@ Data.Ymid=Rangy(2)+(Rangy(1)-Rangy(2))*(Data.Civ3_Y-0.5)/(Npy-1);%temporary coordinate (velocity taken at the point middle from imgae 1 and 2) Data.Uphys=Data.Civ3_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1); - Data.Vphys=Data.Civ3_V_smooth*(Rangy(2)-Rangy(1))/(Npy-1); - [Data.Zphys,Data.Xphys,Data.Yphys,Data.Civ3_E]=shift2z(Data.Xmid,Data.Ymid,Data.Uphys,Data.Vphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytemp) + Data.Vphys=Data.Civ3_V_smooth*(Rangy(1)-Rangy(2))/(Npy-1); + [Data.Zphys,Data.Xphys,Data.Yphys,Data.Error]=shift2z(Data.Xmid,Data.Ymid,Data.Uphys,Data.Vphys,XmlData); %Data.Xphys and Data.Xphys are real coordinate (geometric correction more accurate than xtemp/ytemp) if ~isempty(errormsg) disp_uvmat('ERROR',errormsg,checkrun) @@ -813,7 +807,5 @@ end - - - + %% write result in a netcdf file if requested @@ -829,6 +821,6 @@ else % store only phys data - Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys'}; - Data_light.VarDimName={'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}; + Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys','Error'}; + Data_light.VarDimName={'nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3','nb_vec_3'}; temp=find(Data.Civ3_FF==0); Data_light.Zphys=Data.Zphys(temp); @@ -840,5 +832,5 @@ Data_light.Uphys=Data.Uphys(temp); Data_light.Vphys=Data.Vphys(temp); - + Data_light.Error=Data.Error(temp); if exist('ncfile2','var') errormsg=struct2nc(ncfile2,Data_light); @@ -853,5 +845,5 @@ end - toc + disp(['ellapsed time for the loop ' num2str(toc) ' s']) @@ -912,8 +904,8 @@ %% prepare correlation and search boxes -ibx2=ceil(par_civ.CorrBoxSize(1)/2); -iby2=ceil(par_civ.CorrBoxSize(2)/2); -isx2=ceil(par_civ.SearchBoxSize(1)/2); -isy2=ceil(par_civ.SearchBoxSize(2)/2); +ibx2=floor(par_civ.CorrBoxSize(1)/2); +iby2=floor(par_civ.CorrBoxSize(2)/2); +isx2=floor(par_civ.SearchBoxSize(1)/2); +isy2=floor(par_civ.SearchBoxSize(2)/2); shiftx=round(par_civ.SearchBoxShift(:,1)); shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases @@ -944,23 +936,4 @@ check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma); -% %% prepare images -% if isfield(par_civ,'reverse_pair') -% if par_civ.reverse_pair -% if ischar(par_civ.ImageB) -% temp=par_civ.ImageA; -% par_civ.ImageA=imread(par_civ.ImageB); -% end -% if ischar(temp) -% par_civ.ImageB=imread(temp); -% end -% end -% else -% if ischar(par_civ.ImageA) -% par_civ.ImageA=imread(par_civ.ImageA); -% end -% if ischar(par_civ.ImageB) -% par_civ.ImageB=imread(par_civ.ImageB); -% end -% end par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images par_civ.ImageB=sum(double(par_civ.ImageB),3); @@ -989,6 +962,6 @@ % check_noflux=(par_civ.Mask<100) ;%TODO: to implement check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); - par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area - par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area +% par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area +% par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area end @@ -997,8 +970,7 @@ sum_square=1;% default mesh=1;% default -CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1; -if CheckDecimal - mesh=0.2;%mesh in pixels for subpixel image interpolation - CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; +CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; +if CheckDeformation + mesh=0.25;%mesh in pixels for subpixel image interpolation end % vector=[0 0];%default @@ -1020,13 +992,24 @@ image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A + mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask + mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 - check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; - + check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B - image1_mean=mean(mean(image1_crop)); - image2_mean=mean(mean(image2_crop)); + mask1_crop(check1_y,check1_x)=check_undefined(subrange1_y(check1_y),subrange1_x(check1_x));%extract a mask subimage (correlation box) from image A + mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from imag + sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image + if sizemask > 1/2% eliminate point if more than half of the correlation box is masked + F(ivec)=3; % + else + image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) + image2_crop=image2_crop.*~mask2_crop; + image1_mean=mean(mean(image1_crop))/(1-sizemask); + image2_mean=mean(mean(image2_crop))/(1-sizemask); +% image1_mean=mean(mean(image1_crop)); +% image2_mean=mean(mean(image2_crop)); %threshold on image minimum if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) @@ -1037,22 +1020,20 @@ F(ivec)=3; end - % end + end if F(ivec)~=3 - image1_crop=image1_crop-image1_mean;%substract the mean - image2_crop=image2_crop-image2_mean; - if CheckDecimal - xi=(1:mesh:size(image1_crop,2)); - yi=(1:mesh:size(image1_crop,1))'; - if CheckDeformation + image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts + image2_crop=(image2_crop-image2_mean).*~mask2_crop; + if CheckDeformation + xi=(1:mesh:size(image1_crop,2)); + yi=(1:mesh:size(image1_crop,1))'; [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2)); XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); image1_crop=interp2(image1_crop,XIant,YIant); - else - image1_crop=interp2(image1_crop,xi,yi); - end - xi=(1:mesh:size(image2_crop,2)); - yi=(1:mesh:size(image2_crop,1))'; - image2_crop=interp2(image2_crop,xi,yi); + image1_crop(isnan(image1_crop))=0; + xi=(1:mesh:size(image2_crop,2)); + yi=(1:mesh:size(image2_crop,1))'; + image2_crop=interp2(image2_crop,xi,yi,'*spline'); + image2_crop(isnan(image2_crop))=0; end sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2; @@ -1063,4 +1044,7 @@ %Find the correlation max, at 255 [y,x] = find(result_conv==255,1); + subimage2_crop=image2_crop(y:y+2*iby2/mesh,x:x+2*ibx2/mesh);%subimage of image 2 corresponding to the optimum displacement of first image + sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 + sum_square=sqrt(sum_square);% sqrt of the variance product to normalise correlation if ~isempty(y) && ~isempty(x) try @@ -1286,5 +1270,5 @@ % ymid+ v/2: set of apparent phys y coordinates in the ref plane, image B % XmlData: content of the xml files containing geometric calibration parameters -function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) +function [z,Xphy,Yphy,Error]=shift2z(xmid, ymid, u, v,XmlData) z=0; error=0; @@ -1295,5 +1279,5 @@ R=(Calib_A.R)'; x_a=xmid- u/2; -y_a=ymid- u/2; +y_a=ymid- v/2; z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3); Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a; @@ -1311,6 +1295,10 @@ %% second image +%loading shift angle + Calib_B=XmlData{2}.GeometryCalib; R=(Calib_B.R)'; + + x_b=xmid+ u/2; y_b=ymid+ v/2; @@ -1330,5 +1318,11 @@ %% result Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya); -error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den; +mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; +mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; +mtz=(XmlData{1}.GeometryCalib.Tx_Ty_Tz(1,3)+XmlData{2}.GeometryCalib.Tx_Ty_Tz(1,3))/2; + +Error=(sqrt(mfx^2+mfy^2)/(2*sqrt(2)*mtz)).*(((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den); + +z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den; xnew(1,:)=Dxa.*z+x_a; @@ -1336,10 +1330,9 @@ ynew(1,:)=Dya.*z+y_a; ynew(2,:)=Dyb.*z+y_b; - -Xphy=mean(xnew,1); %on moyenne les 2 valeurs +Xphy=mean(xnew,1); Yphy=mean(ynew,1); -z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den; - - - - + + + + + Index: /trunk/src/tutorial_jet.m =================================================================== --- /trunk/src/tutorial_jet.m (revision 878) +++ /trunk/src/tutorial_jet.m (revision 878) @@ -0,0 +1,18 @@ +DataFolder=''; %TODO: put the actual path to start browser +DataFolder='C:\Users\sommeria\Documents\MATLAB_WORK\UVMAT_DEMO_SOURCES\UVMAT_DEMO03_PIVchallenge_2005C'; +fileinput=uigetfile_uvmat('pick an input file',DataFolder); +Data=nc2struct(fileinput); +figure +npy=numel(Data.coord_y); +npx=numel(Data.coord_x); +X=ones(npy,1)*Data.coord_x'; +Sum=sum(Data.VMean,2); +XMean=sum(X.*Data.VMean,2)./Sum; +XMean=XMean*ones(1,npx); +Xrms=sum(((X-XMean).*(X-XMean).*Data.VMean),2)./Sum; +Xrm=sqrt(Xrms); +plot(Data.coord_y(3:end-1),Xrms(3:end-1)) +%0.024 y+0.14 +% spread= sqrt(2*log(2))*0.024=0.028 +%y0=-5.8; +