Changeset 878


Ignore:
Timestamp:
Feb 23, 2015, 10:17:04 PM (10 years ago)
Author:
sommeria
Message:

stereo_civ and civ2vel_3C updated, spline introduced for civ_series:deformation

Location:
trunk/src
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/civ2vel_3C.m

    r864 r878  
    169169%% MAIN LOOP ON FIELDS
    170170warning off
    171 if NbField<2 && length(filecell{:,1})>2
     171if NbField<2
    172172    disp_uvmat('ERROR','you need at least 2 images to compute the mean position for the stereo.',checkrun)
    173173return
    … …  
    285285end
    286286    ZI=mean(ZItemp,3); %mean between two the two time step
     287    Vtest=ZItemp(:,:,2)-ZItemp(:,:,1);
    287288   
    288289    [Xa,Ya]=px_XYZ(XmlData{1}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view a
    289290    [Xb,Yb]=px_XYZ(XmlData{2}.GeometryCalib,XI,YI,ZI);% set of image coordinates on view b
    290291   
    291     %trouver z
    292     % trouver les coordonnées px sur chaque image.
    293     %A=
     292   
    294293    for iview=1:2
    295294        %% reading input file(s)
    … …  
    322321    Ua=griddata(X1,Y1,U1,Xa,Ya);
    323322    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
    327328    temp=find(Data{2}.FF==0);
    328329    X2=Data{2}.X(temp);
    … …  
    332333    Ub=griddata(X2,Y2,U2,Xb,Yb);
    333334    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
    335341    S=ones(size(XI,1),size(XI,2),3);
    336342    D=ones(size(XI,1),size(XI,2),3,3);
     343
    337344    S(:,:,1)=A(:,:,1,1).*Ua+A(:,:,2,1).*Va+B(:,:,1,1).*Ub+B(:,:,2,1).*Vb;
    338345    S(:,:,2)=A(:,:,1,2).*Ua+A(:,:,2,2).*Va+B(:,:,1,2).*Ub+B(:,:,2,2).*Vb;
    … …  
    349356    for indj=1:size(XI,1)
    350357        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...
    352359            U(indj,indi)=dxyz(1);
    353360            V(indj,indi)=dxyz(2);
    … …  
    360367    Error(:,:,3)=B(:,:,1,1).*U+B(:,:,1,2).*V+B(:,:,1,3).*W-Ub;
    361368    Error(:,:,4)=B(:,:,2,1).*U+B(:,:,2,2).*V+B(:,:,2,3).*W-Vb;
     369   
     370   
     371
    362372   
    363373    %% generating the name of the merged field
    … …  
    396406    MergeData.V=V/Dt;
    397407    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));
    399412    errormsg=struct2nc(OutputFile,MergeData);%save result file
    400413    if isempty(errormsg)
    … …  
    406419
    407420
    408 function A=get_coeff(Calib,X,Y,x,y,z)
     421function [A]=get_coeff(Calib,X,Y,x,y,z) % compute A~ coefficients
    409422R=(Calib.R)';%rotation matrix
    410423T_z=Calib.Tx_Ty_Tz(3);
    411424T=R(7)*x+R(8)*y+R(9)*z+T_z;
     425
    412426A(:,:,1,1)=(R(1)-R(7)*X)./T;
    413427A(:,:,1,2)=(R(2)-R(8)*X)./T;
    … …  
    417431A(:,:,2,3)=(R(6)-R(9)*Y)./T;
    418432
    419 
    420 
    421 
    422 function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData)
     433function [U,V,X,Y]=Ud2U(Calib,Xd,Yd,Ud,Vd) % convert Xd to X  and Ud to U
     434
     435X1d=Xd-Ud/2;
     436X2d=Xd+Ud/2;
     437Y1d=Yd-Vd/2;
     438Y2d=Yd+Vd/2;
     439
     440X1=(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);
     441X2=(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);
     442Y1=(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);
     443Y2=(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
     445U=X2-X1;
     446V=Y2-Y1;
     447X=X1+U/2;
     448Y=Y1+V/2;
     449
     450
     451
     452function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData) % get H from stereo data
    423453z=0;
    424454error=0;
    … …  
    429459R=(Calib_A.R)';
    430460x_a=xmid- u/2;
    431 y_a=ymid- u/2;
     461y_a=ymid- v/2;
    432462z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3);
    433463Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a;
    … …  
    445475
    446476%% second image
     477%loading shift angle
     478
    447479Calib_B=XmlData{2}.GeometryCalib;
    448480R=(Calib_B.R)';
     481
     482
    449483x_b=xmid+ u/2;
    450484y_b=ymid+ v/2;
    … …  
    464498%% result
    465499Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
    466 error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
     500error=((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);
     506z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den;
    467507
    468508xnew(1,:)=Dxa.*z+x_a;
    … …  
    470510ynew(1,:)=Dya.*z+y_a;
    471511ynew(2,:)=Dyb.*z+y_b;
    472 
    473 Xphy=mean(xnew,1); %on moyenne les 2 valeurs
    474 Yphy=mean(ynew,1);
    475 z=((Dxb-Dxa).*u-(Dyb-Dya).*v)./Den;
    476 
    477 
     512Xphy=mean(xnew,1);
     513Yphy=mean(ynew,1);
     514
     515
     516
     517

  • trunk/src/series/civ_series.m

    r876 r878  
    324324%%%%% MAIN LOOP %%%%%%
    325325maskname='';% initiate the mask name
     326tic;
    326327for ifield=1:NbField
    327328    if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
    … …  
    852853    end
    853854end
     855disp(['ellapsed time ' num2str(toc) ' s'])
    854856
    855857% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
    … …  
    10241026                XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    10251027                YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
    1026                         image1_crop_1=image1_crop;
    10271028                image1_crop=interp2(image1_crop,XIant,YIant);
    1028                         image1_crop_2=image1_crop;
    10291029                image1_crop(isnan(image1_crop))=0;
    1030                   image1_crop_3=image1_crop;
    10311030                xi=(1:mesh:size(image2_crop,2));
    10321031                yi=(1:mesh:size(image2_crop,1))';
    1033                 image2_crop=interp2(image2_crop,xi,yi);
     1032                image2_crop=interp2(image2_crop,xi,yi,'*spline');
    10341033                image2_crop(isnan(image2_crop))=0;
    1035                 par_civ.CorrSmooth=3;
     1034                %par_civ.CorrSmooth=3;%%%%%%%%%%%%%%%%%%%
    10361035                %%
    10371036            end

  • trunk/src/series/stereo_civ.m

    r876 r878  
    1 %'stereo_civ': determination of topography by image correlation of two stereo views
     1%'civ_series': PIV function activated by the general GUI series
    22% --- call the sub-functions:
    33%   civ: PIV function itself
    … …  
    282282       
    283283    %% Civ1
    284    
     284
     285   
    285286    % if Civ1 computation is requested
    286287    if isfield (Param.ActionInput,'Civ1')
    … …  
    297298            end
    298299        end
     300       
     301
     302       
    299303        [A,Rangx,Rangy]=phys_ima(A,XmlData,1);
    300304        [Npy,Npx]=size(A{1});
    … …  
    516520            mask=imread(par_civ2.Mask);
    517521        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);
    522524        par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
    523525        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
    524526        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);
    529531        end
    530532        % calculate velocity data (y and v in indices, reverse to y component)
    … …  
    617619        Data.Civ2_FF(ind_good)=FFres;
    618620        Data.CivStage=Data.CivStage+1;
    619        
    620            
    621  
    622        
    623621    end
    624622   
    … …  
    659657            [GridX,GridY]=meshgrid(minix:par_civ3.Dx:maxix,miniy:par_civ3.Dy:maxiy);
    660658            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);       
    664660        end
    665661        Shiftx=zeros(size(par_civ3.Grid,1),1);% shift expected from civ2 data
    … …  
    696692            mask=imread(par_civ3.Mask);
    697693        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);
    702696        par_civ3.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
    703697        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
    704698        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);
    709703        end
    710704        % calculate velocity data (y and v in indices, reverse to y component)
    … …  
    724718       
    725719        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 record
    727         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'}];
    728722        Data.VarAttribute{nbvar+1}.Role='coord_x';
    729723        Data.VarAttribute{nbvar+2}.Role='coord_y';
    … …  
    776770        Data.Patch3_SubDomainSize=Param.ActionInput.Patch3.SubDomainSize;
    777771        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'}];
    779773        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'}];
    781775       
    782776        Data.VarAttribute{nbvar+1}.Role='vector_x';
    … …  
    805799        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)
    806800        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)
    809803        if ~isempty(errormsg)
    810804            disp_uvmat('ERROR',errormsg,checkrun)
    … …  
    813807       
    814808    end
    815    
    816    
    817    
     809     
    818810   
    819811    %% write result in a netcdf file if requested
    … …  
    829821    else
    830822       % 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'};
    833825        temp=find(Data.Civ3_FF==0);
    834826        Data_light.Zphys=Data.Zphys(temp);
    … …  
    840832        Data_light.Uphys=Data.Uphys(temp);
    841833        Data_light.Vphys=Data.Vphys(temp);
    842        
     834        Data_light.Error=Data.Error(temp);
    843835       if exist('ncfile2','var')
    844836            errormsg=struct2nc(ncfile2,Data_light);
    … …  
    853845   end
    854846
    855   toc
     847  disp(['ellapsed time for the loop ' num2str(toc) ' s'])
    856848
    857849
    … …  
    912904
    913905%% 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);
     906ibx2=floor(par_civ.CorrBoxSize(1)/2);
     907iby2=floor(par_civ.CorrBoxSize(2)/2);
     908isx2=floor(par_civ.SearchBoxSize(1)/2);
     909isy2=floor(par_civ.SearchBoxSize(2)/2);
    918910shiftx=round(par_civ.SearchBoxShift(:,1));
    919911shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
    … …  
    944936check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
    945937
    946 % %% prepare images
    947 % if isfield(par_civ,'reverse_pair')
    948 %     if par_civ.reverse_pair
    949 %         if ischar(par_civ.ImageB)
    950 %             temp=par_civ.ImageA;
    951 %             par_civ.ImageA=imread(par_civ.ImageB);
    952 %         end
    953 %         if ischar(temp)
    954 %             par_civ.ImageB=imread(temp);
    955 %         end
    956 %     end
    957 % else
    958 %     if ischar(par_civ.ImageA)
    959 %         par_civ.ImageA=imread(par_civ.ImageA);
    960 %     end
    961 %     if ischar(par_civ.ImageB)
    962 %         par_civ.ImageB=imread(par_civ.ImageB);
    963 %     end
    964 % end
    965938par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
    966939par_civ.ImageB=sum(double(par_civ.ImageB),3);
    … …  
    989962  %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
    990963    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  area
    992     par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
     964%     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
    993966end
    994967
    … …  
    997970sum_square=1;% default
    998971mesh=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;
     972CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
     973if CheckDeformation
     974    mesh=0.25;%mesh in pixels for subpixel image interpolation
    1003975end
    1004976% vector=[0 0];%default
    … …  
    1020992    image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    1021993    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
    1022996    check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    1023997    check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    1024998    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;   
    10271000    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
    10281001    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));
    10311014    %threshold on image minimum
    10321015    if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    … …  
    10371020        F(ivec)=3;
    10381021    end
    1039     %         end
     1022    end
    10401023    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))';
    10471029                [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
    10481030                XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    10491031                YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
    10501032                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;
    10571038        end
    10581039        sum_square=(sum(sum(image1_crop.*image1_crop)));%+sum(sum(image2_crop.*image2_crop)))/2;
    … …  
    10631044        %Find the correlation max, at 255
    10641045        [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
    10651049        if ~isempty(y) && ~isempty(x)
    10661050            try
    … …  
    12861270% ymid+ v/2: set of apparent phys y coordinates in the ref plane, image B
    12871271% XmlData: content of the xml files containing geometric calibration parameters
    1288 function [z,Xphy,Yphy,error]=shift2z(xmid, ymid, u, v,XmlData)
     1272function [z,Xphy,Yphy,Error]=shift2z(xmid, ymid, u, v,XmlData)
    12891273z=0;
    12901274error=0;
    … …  
    12951279R=(Calib_A.R)';
    12961280x_a=xmid- u/2;
    1297 y_a=ymid- u/2;
     1281y_a=ymid- v/2;
    12981282z_a=R(7)*x_a+R(8)*y_a+Calib_A.Tx_Ty_Tz(1,3);
    12991283Xa=(R(1)*x_a+R(2)*y_a+Calib_A.Tx_Ty_Tz(1,1))./z_a;
    … …  
    13111295
    13121296%% second image
     1297%loading shift angle
     1298
    13131299Calib_B=XmlData{2}.GeometryCalib;
    13141300R=(Calib_B.R)';
     1301
     1302
    13151303x_b=xmid+ u/2;
    13161304y_b=ymid+ v/2;
    … …  
    13301318%% result
    13311319Den=(Dxb-Dxa).*(Dxb-Dxa)+(Dyb-Dya).*(Dyb-Dya);
    1332 error=((Dyb-Dya).*u-(Dxb-Dxa).*v)./Den;
     1320mfx=(XmlData{1}.GeometryCalib.fx_fy(1)+XmlData{2}.GeometryCalib.fx_fy(1))/2;
     1321mfy=(XmlData{1}.GeometryCalib.fx_fy(2)+XmlData{2}.GeometryCalib.fx_fy(2))/2;
     1322mtz=(XmlData{1}.GeometryCalib.Tx_Ty_Tz(1,3)+XmlData{2}.GeometryCalib.Tx_Ty_Tz(1,3))/2;
     1323
     1324Error=(sqrt(mfx^2+mfy^2)/(2*sqrt(2)*mtz)).*(((Dyb-Dya).*(-u)-(Dxb-Dxa).*(-v))./Den);
     1325
     1326z=((Dxb-Dxa).*(-u)+(Dyb-Dya).*(-v))./Den;
    13331327
    13341328xnew(1,:)=Dxa.*z+x_a;
    … …  
    13361330ynew(1,:)=Dya.*z+y_a;
    13371331ynew(2,:)=Dyb.*z+y_b;
    1338 
    1339 Xphy=mean(xnew,1); %on moyenne les 2 valeurs
     1332Xphy=mean(xnew,1);
    13401333Yphy=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.