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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.