Changeset 878 for trunk/src/series
- Timestamp:
- Feb 23, 2015, 10:17:04 PM (10 years ago)
- Location:
- trunk/src/series
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/series/civ2vel_3C.m
r864 r878 169 169 %% MAIN LOOP ON FIELDS 170 170 warning off 171 if NbField<2 && length(filecell{:,1})>2171 if NbField<2 172 172 disp_uvmat('ERROR','you need at least 2 images to compute the mean position for the stereo.',checkrun) 173 173 return … … 285 285 end 286 286 ZI=mean(ZItemp,3); %mean between two the two time step 287 Vtest=ZItemp(:,:,2)-ZItemp(:,:,1); 287 288 288 289 [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view a 289 290 [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b 290 291 291 %trouver z 292 % trouver les coordonnées px sur chaque image. 293 %A= 292 294 293 for iview=1:2 295 294 %% reading input file(s) … … 322 321 Ua=griddata(X1,Y1,U1,Xa,Ya); 323 322 Va=griddata(X1,Y1,V1,Xa,Ya); 324 A=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,YI,YI,ZI); 325 326 323 324 [Ua,Va,Xa,Ya]=Ud2U(XmlData{1}.GeometryCalib,Xa,Ya,Ua,Va); % convert Xd data to X 325 [A]=get_coeff(XmlData{1}.GeometryCalib,Xa,Ya,XI,YI,ZI); %get coef A~ 326 327 %remove wrong vector 327 328 temp=find(Data{2}.FF==0); 328 329 X2=Data{2}.X(temp); … … 332 333 Ub=griddata(X2,Y2,U2,Xb,Yb); 333 334 Vb=griddata(X2,Y2,V2,Xb,Yb); 334 B=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,YI,YI,ZI); 335 336 [Ub,Vb,Xb,Yb]=Ud2U(XmlData{2}.GeometryCalib,Xb,Yb,Ub,Vb); % convert Xd data to X 337 [B]=get_coeff(XmlData{2}.GeometryCalib,Xb,Yb,XI,YI,ZI); %get coef B~ 338 339 340 % System to solve 335 341 S=ones(size(XI,1),size(XI,2),3); 336 342 D=ones(size(XI,1),size(XI,2),3,3); 343 337 344 S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb; 338 345 S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb; … … 349 356 for indj=1:size(XI,1) 350 357 for indi=1:size(XI,2) 351 dxyz= squeeze(S(indj,indi,:))\squeeze(D(indj,indi,:,:));358 dxyz=(squeeze(D(indj,indi,:,:))*1000)\(squeeze(S(indj,indi,:))*1000); % solving... 352 359 U(indj,indi)=dxyz(1); 353 360 V(indj,indi)=dxyz(2); … … 360 367 Error(:,:,3)=B(:,:,1,1).*U+B(:,:,1,2).*V+B(:,:,1,3).*W-Ub; 361 368 Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb; 369 370 371 362 372 363 373 %% generating the name of the merged field … … 396 406 MergeData.V=V/Dt; 397 407 MergeData.W=W/Dt; 398 MergeData.Error=sqrt(sum(Error.*Error,3)); 408 409 mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; 410 mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; 411 MergeData.Error=(sqrt(mfx^2+mfy^2)/4).*sqrt(sum(Error.*Error,3)); 399 412 errormsg=struct2nc(OutputFile,MergeData);%save result file 400 413 if isempty(errormsg) … … 406 419 407 420 408 function A=get_coeff(Calib,X,Y,x,y,z)421 function [A]=get_coeff(Calib,X,Y,x,y,z) % compute A~ coefficients 409 422 R=(Calib.R)';%rotation matrix 410 423 T_z=Calib.Tx_Ty_Tz(3); 411 424 T=R(7)*x+R(8)*y+R(9)*z+T_z; 425 412 426 A(:,:,1,1)=(R(1)-R(7)*X)./T; 413 427 A(:,:,1,2)=(R(2)-R(8)*X)./T; … … 417 431 A(:,:,2,3)=(R(6)-R(9)*Y)./T; 418 432 419 420 421 422 function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) 433 function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X and Ud to U 434 435 X1d=Xd-Ud/2; 436 X2d=Xd+Ud/2; 437 Y1d=Yd-Vd/2; 438 Y2d=Yd+Vd/2; 439 440 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); 441 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); 442 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); 443 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); 444 445 U=X2-X1; 446 V=Y2-Y1; 447 X=X1+U/2; 448 Y=Y1+V/2; 449 450 451 452 function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) % get H from stereo data 423 453 z=0; 424 454 error=0; … … 429 459 R=(Calib_A.R)'; 430 460 x_a=xmid- u/2; 431 y_a=ymid- u/2;461 y_a=ymid- v/2; 432 462 z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3); 433 463 Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a; … … 445 475 446 476 %% second image 477 %loading shift angle 478 447 479 Calib_B=XmlData{2}.GeometryCalib; 448 480 R=(Calib_B.R)'; 481 482 449 483 x_b=xmid+ u/2; 450 484 y_b=ymid+ v/2; … … 464 498 %% result 465 499 Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya); 466 error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den; 500 error=((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den; 501 % ex=-error.*(Dyb-Dya); 502 % ey=-error.*(Dxb-Dxa); 503 504 % z1=-u./(Dxb-Dxa); 505 % z2=-v./(Dyb-Dya); 506 z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den; 467 507 468 508 xnew(1,:)=Dxa.*z+x_a; … … 470 510 ynew(1,:)=Dya.*z+y_a; 471 511 ynew(2,:)=Dyb.*z+y_b; 472 473 Xphy=mean(xnew,1); %on moyenne les 2 valeurs474 Yphy=mean(ynew,1); 475 z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den; 476 477 512 Xphy=mean(xnew,1); 513 Yphy=mean(ynew,1); 514 515 516 517 -
trunk/src/series/civ_series.m
r876 r878 324 324 %%%%% MAIN LOOP %%%%%% 325 325 maskname='';% initiate the mask name 326 tic; 326 327 for ifield=1:NbField 327 328 if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series (checkrun=1) … … 852 853 end 853 854 end 855 disp(['ellapsed time ' num2str(toc) ' s']) 854 856 855 857 % 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/ … … 1024 1026 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 1025 1027 YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); 1026 image1_crop_1=image1_crop;1027 1028 image1_crop=interp2(image1_crop,XIant,YIant); 1028 image1_crop_2=image1_crop;1029 1029 image1_crop(isnan(image1_crop))=0; 1030 image1_crop_3=image1_crop;1031 1030 xi=(1:mesh:size(image2_crop,2)); 1032 1031 yi=(1:mesh:size(image2_crop,1))'; 1033 image2_crop=interp2(image2_crop,xi,yi );1032 image2_crop=interp2(image2_crop,xi,yi,'*spline'); 1034 1033 image2_crop(isnan(image2_crop))=0; 1035 par_civ.CorrSmooth=3;1034 %par_civ.CorrSmooth=3;%%%%%%%%%%%%%%%%%%% 1036 1035 %% 1037 1036 end -
trunk/src/series/stereo_civ.m
r876 r878 1 %' stereo_civ': determination of topography by image correlation of two stereo views1 %'civ_series': PIV function activated by the general GUI series 2 2 % --- call the sub-functions: 3 3 % civ: PIV function itself … … 282 282 283 283 %% Civ1 284 284 285 285 286 % if Civ1 computation is requested 286 287 if isfield (Param.ActionInput,'Civ1') … … 297 298 end 298 299 end 300 301 302 299 303 [A,Rangx,Rangy]=phys_ima(A,XmlData,1); 300 304 [Npy,Npx]=size(A{1}); … … 516 520 mask=imread(par_civ2.Mask); 517 521 end 518 ibx2=ceil(par_civ2.CorrBoxSize(1)/2); 519 iby2=ceil(par_civ2.CorrBoxSize(2)/2); 520 % par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess 521 % par_civ2.SearchBoxSize(2)=2*iby2+9; 522 % ibx2=ceil(par_civ2.CorrBoxSize(1)/2); 523 % iby2=ceil(par_civ2.CorrBoxSize(2)/2); 522 524 par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; 523 525 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 524 526 if par_civ2.CheckDeformation 525 par_civ2.DUDX=DUDX ./nbval;526 par_civ2.DUDY=DUDY ./nbval;527 par_civ2.DVDX=DVDX ./nbval;528 par_civ2.DVDY=DVDY ./nbval;527 par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1); 528 par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1); 529 par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1); 530 par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); 529 531 end 530 532 % calculate velocity data (y and v in indices, reverse to y component) … … 617 619 Data.Civ2_FF(ind_good)=FFres; 618 620 Data.CivStage=Data.CivStage+1; 619 620 621 622 623 621 end 624 622 … … 659 657 [GridX,GridY]=meshgrid(minix:par_civ3.Dx:maxix,miniy:par_civ3.Dy:maxiy); 660 658 par_civ3.Grid(:,1)=reshape(GridX,[],1); 661 par_civ3.Grid(:,2)=reshape(GridY,[],1); 662 663 659 par_civ3.Grid(:,2)=reshape(GridY,[],1); 664 660 end 665 661 Shiftx=zeros(size(par_civ3.Grid,1),1);% shift expected from civ2 data … … 696 692 mask=imread(par_civ3.Mask); 697 693 end 698 ibx2=ceil(par_civ3.CorrBoxSize(1)/2); 699 iby2=ceil(par_civ3.CorrBoxSize(2)/2); 700 % par_civ3.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess 701 % par_civ3.SearchBoxSize(2)=2*iby2+9; 694 % ibx2=ceil(par_civ3.CorrBoxSize(1)/2); 695 % iby2=ceil(par_civ3.CorrBoxSize(2)/2); 702 696 par_civ3.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)]; 703 697 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 704 698 if par_civ3.CheckDeformation 705 par_civ3.DUDX=DUDX ./nbval;706 par_civ3.DUDY=DUDY ./nbval;707 par_civ3.DVDX=DVDX ./nbval;708 par_civ3.DVDY=DVDY ./nbval;699 par_civ3.DUDX=DUDX(nbval>=1)./nbval(nbval>=1); 700 par_civ3.DUDY=DUDY(nbval>=1)./nbval(nbval>=1); 701 par_civ3.DVDX=DVDX(nbval>=1)./nbval(nbval>=1); 702 par_civ3.DVDY=DVDY(nbval>=1)./nbval(nbval>=1); 709 703 end 710 704 % calculate velocity data (y and v in indices, reverse to y component) … … 724 718 725 719 nbvar=numel(Data.ListVarName); 726 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 record727 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'}];720 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 721 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'}]; 728 722 Data.VarAttribute{nbvar+1}.Role='coord_x'; 729 723 Data.VarAttribute{nbvar+2}.Role='coord_y'; … … 776 770 Data.Patch3_SubDomainSize=Param.ActionInput.Patch3.SubDomainSize; 777 771 nbvar=length(Data.ListVarName); 778 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' }];772 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'}]; 779 773 Data.VarDimName=[Data.VarDimName {'nb_vec_3','nb_vec_3',{'nb_coord','nb_bounds','nb_subdomain_3'},{'nb_subdomain_3'},... 780 {'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' }];774 {'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'}]; 781 775 782 776 Data.VarAttribute{nbvar+1}.Role='vector_x'; … … 805 799 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) 806 800 Data.Uphys=Data.Civ3_U_smooth*(Rangx(2)-Rangx(1))/(Npx-1); 807 Data.Vphys=Data.Civ3_V_smooth*(Rangy( 2)-Rangy(1))/(Npy-1);808 [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)801 Data.Vphys=Data.Civ3_V_smooth*(Rangy(1)-Rangy(2))/(Npy-1); 802 [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) 809 803 if ~isempty(errormsg) 810 804 disp_uvmat('ERROR',errormsg,checkrun) … … 813 807 814 808 end 815 816 817 809 818 810 819 811 %% write result in a netcdf file if requested … … 829 821 else 830 822 % store only phys data 831 Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys' };832 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' };823 Data_light.ListVarName={'Xphys','Yphys','Zphys','Civ3_C','Xmid','Ymid','Uphys','Vphys','Error'}; 824 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'}; 833 825 temp=find(Data.Civ3_FF==0); 834 826 Data_light.Zphys=Data.Zphys(temp); … … 840 832 Data_light.Uphys=Data.Uphys(temp); 841 833 Data_light.Vphys=Data.Vphys(temp); 842 834 Data_light.Error=Data.Error(temp); 843 835 if exist('ncfile2','var') 844 836 errormsg=struct2nc(ncfile2,Data_light); … … 853 845 end 854 846 855 toc847 disp(['ellapsed time for the loop ' num2str(toc) ' s']) 856 848 857 849 … … 912 904 913 905 %% prepare correlation and search boxes 914 ibx2= ceil(par_civ.CorrBoxSize(1)/2);915 iby2= ceil(par_civ.CorrBoxSize(2)/2);916 isx2= ceil(par_civ.SearchBoxSize(1)/2);917 isy2= ceil(par_civ.SearchBoxSize(2)/2);906 ibx2=floor(par_civ.CorrBoxSize(1)/2); 907 iby2=floor(par_civ.CorrBoxSize(2)/2); 908 isx2=floor(par_civ.SearchBoxSize(1)/2); 909 isy2=floor(par_civ.SearchBoxSize(2)/2); 918 910 shiftx=round(par_civ.SearchBoxShift(:,1)); 919 911 shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases … … 944 936 check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma); 945 937 946 % %% prepare images947 % if isfield(par_civ,'reverse_pair')948 % if par_civ.reverse_pair949 % if ischar(par_civ.ImageB)950 % temp=par_civ.ImageA;951 % par_civ.ImageA=imread(par_civ.ImageB);952 % end953 % if ischar(temp)954 % par_civ.ImageB=imread(temp);955 % end956 % end957 % else958 % if ischar(par_civ.ImageA)959 % par_civ.ImageA=imread(par_civ.ImageA);960 % end961 % if ischar(par_civ.ImageB)962 % par_civ.ImageB=imread(par_civ.ImageB);963 % end964 % end965 938 par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images 966 939 par_civ.ImageB=sum(double(par_civ.ImageB),3); … … 989 962 % check_noflux=(par_civ.Mask<100) ;%TODO: to implement 990 963 check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 ); 991 par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area992 par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area964 % par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined area 965 % par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined area 993 966 end 994 967 … … 997 970 sum_square=1;% default 998 971 mesh=1;% default 999 CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1; 1000 if CheckDecimal 1001 mesh=0.2;%mesh in pixels for subpixel image interpolation 1002 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; 972 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; 973 if CheckDeformation 974 mesh=0.25;%mesh in pixels for subpixel image interpolation 1003 975 end 1004 976 % vector=[0 0];%default … … 1020 992 image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A 1021 993 image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A 994 mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask 995 mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask 1022 996 check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1 1023 997 check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight; 1024 998 check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2 1025 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 1026 999 check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight; 1027 1000 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 1028 1001 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 1029 image1_mean=mean(mean(image1_crop)); 1030 image2_mean=mean(mean(image2_crop)); 1002 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 1003 mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from imag 1004 sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image 1005 if sizemask > 1/2% eliminate point if more than half of the correlation box is masked 1006 F(ivec)=3; % 1007 else 1008 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1) 1009 image2_crop=image2_crop.*~mask2_crop; 1010 image1_mean=mean(mean(image1_crop))/(1-sizemask); 1011 image2_mean=mean(mean(image2_crop))/(1-sizemask); 1012 % image1_mean=mean(mean(image1_crop)); 1013 % image2_mean=mean(mean(image2_crop)); 1031 1014 %threshold on image minimum 1032 1015 if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma) … … 1037 1020 F(ivec)=3; 1038 1021 end 1039 %end1022 end 1040 1023 if F(ivec)~=3 1041 image1_crop=image1_crop-image1_mean;%substract the mean 1042 image2_crop=image2_crop-image2_mean; 1043 if CheckDecimal 1044 xi=(1:mesh:size(image1_crop,2)); 1045 yi=(1:mesh:size(image1_crop,1))'; 1046 if CheckDeformation 1024 image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts 1025 image2_crop=(image2_crop-image2_mean).*~mask2_crop; 1026 if CheckDeformation 1027 xi=(1:mesh:size(image1_crop,2)); 1028 yi=(1:mesh:size(image1_crop,1))'; 1047 1029 [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2)); 1048 1030 XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 1049 1031 YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2); 1050 1032 image1_crop=interp2(image1_crop,XIant,YIant); 1051 else 1052 image1_crop=interp2(image1_crop,xi,yi); 1053 end 1054 xi=(1:mesh:size(image2_crop,2)); 1055 yi=(1:mesh:size(image2_crop,1))'; 1056 image2_crop=interp2(image2_crop,xi,yi); 1033 image1_crop(isnan(image1_crop))=0; 1034 xi=(1:mesh:size(image2_crop,2)); 1035 yi=(1:mesh:size(image2_crop,1))'; 1036 image2_crop=interp2(image2_crop,xi,yi,'*spline'); 1037 image2_crop(isnan(image2_crop))=0; 1057 1038 end 1058 1039 sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2; … … 1063 1044 %Find the correlation max, at 255 1064 1045 [y,x] = find(result_conv==255,1); 1046 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 1047 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2 1048 sum_square=sqrt(sum_square);% sqrt of the variance product to normalise correlation 1065 1049 if ~isempty(y) && ~isempty(x) 1066 1050 try … … 1286 1270 % ymid+ v/2: set of apparent phys y coordinates in the ref plane, image B 1287 1271 % XmlData: content of the xml files containing geometric calibration parameters 1288 function [z,Xphy,Yphy, error]=shift2z(xmid, ymid, u, v,XmlData)1272 function [z,Xphy,Yphy,Error]=shift2z(xmid, ymid, u, v,XmlData) 1289 1273 z=0; 1290 1274 error=0; … … 1295 1279 R=(Calib_A.R)'; 1296 1280 x_a=xmid- u/2; 1297 y_a=ymid- u/2;1281 y_a=ymid- v/2; 1298 1282 z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3); 1299 1283 Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a; … … 1311 1295 1312 1296 %% second image 1297 %loading shift angle 1298 1313 1299 Calib_B=XmlData{2}.GeometryCalib; 1314 1300 R=(Calib_B.R)'; 1301 1302 1315 1303 x_b=xmid+ u/2; 1316 1304 y_b=ymid+ v/2; … … 1330 1318 %% result 1331 1319 Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya); 1332 error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den; 1320 mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2; 1321 mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2; 1322 mtz=(XmlData{1}.GeometryCalib.Tx_Ty_Tz(1,3)+XmlData{2}.GeometryCalib.Tx_Ty_Tz(1,3))/2; 1323 1324 Error=(sqrt(mfx^2+mfy^2)/(2*sqrt(2)*mtz)).*(((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den); 1325 1326 z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den; 1333 1327 1334 1328 xnew(1,:)=Dxa.*z+x_a; … … 1336 1330 ynew(1,:)=Dya.*z+y_a; 1337 1331 ynew(2,:)=Dyb.*z+y_b; 1338 1339 Xphy=mean(xnew,1); %on moyenne les 2 valeurs 1332 Xphy=mean(xnew,1); 1340 1333 Yphy=mean(ynew,1); 1341 z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den; 1342 1343 1344 1345 1334 1335 1336 1337 1338
Note: See TracChangeset
for help on using the changeset viewer.