Ignore:
Timestamp:
Jul 29, 2024, 9:43:17 AM (2 months ago)
Author:
sommeria
Message:

civ3D updated

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/civ_3D.m

    r1163 r1164  
    1 %'civ_3D': 3D PIV from image scan in a volume 
     1%'civ_3D': 3D PIV from image scan in a volume
    22
    33%------------------------------------------------------------------------
    … …  
    1212% Param: Matlab structure of input  parameters
    1313%     Param contains info of the GUI series using the fct read_GUI.
    14 %     Param.Action.RUN = 0 (to set the status of the GUI series) or =1 to RUN the computation 
     14%     Param.Action.RUN = 0 (to set the status of the GUI series) or =1 to RUN the computation
    1515%     Param.InputTable: sets the input file(s)
    1616%           if absent, the fct looks for input data in Param.ActionInput     (test mode)
    … …  
    4949    path_series=fileparts(which('series'));
    5050    addpath(fullfile(path_series,'series'))
    51 
    52    Data=civ_input(Param)
    53 
     51   
     52    Data=civ_input(Param)
     53   
    5454    Data.Program=mfilename;%gives the name of the current function
    5555    Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    … …  
    8383        end
    8484    end
    85 
    86 
     85   
     86   
    8787    return
    8888end
    … …  
    114114if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')
    115115    OutputDir=[Param.OutputSubDir Param.OutputDirExt];
    116      OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
     116    OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
    117117else
    118118    disp_uvmat('ERROR','no output folder defined',checkrun)
    … …  
    129129end
    130130
    131 [~,i1_series,~,j1_series,~]=get_file_series(Param);
    132 iview_A=0;% series index (iview) for the first image series
    133 iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
    134 if Param.ActionInput.CheckCiv1
    135     iview_A=1;% usual PIV, the image series is on the first line of the table
    136 elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line
    137     iview_A=2;% the second line is used for the input images of Civ2
    138 end
    139 if iview_A~=0
    140     RootPath_A=Param.InputTable{iview_A,1};
    141     RootFile_A=Param.InputTable{iview_A,3};
    142     SubDir_A=Param.InputTable{iview_A,2};
    143     NomType_A=Param.InputTable{iview_A,4};
    144     FileExt_A=Param.InputTable{iview_A,5};
    145     if iview_B==0
    146         iview_B=iview_A;% the second image series is the same as the first
    147     end
    148     RootPath_B=Param.InputTable{iview_B,1};
    149     RootFile_B=Param.InputTable{iview_B,3};
    150     SubDir_B=Param.InputTable{iview_B,2};
    151     NomType_B=Param.InputTable{iview_B,4};
    152     FileExt_B=Param.InputTable{iview_B,5};
    153 end
     131[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     132
     133if size(Param.InputTable,1)==2% netcdf file as first input line
     134    RootPath_nc=Param.InputTable{1,1};
     135    RootFile_nc=Param.InputTable{1,3};
     136    SubDir_nc=Param.InputTable{1,2};
     137    NomType_nc=Param.InputTable{1,4};
     138    FileExt_nc=Param.InputTable{1,5};
     139    iview_A=2;%second line used for image input
     140    iview_B=2;
     141else
     142    iview_A=1;%second line used for image input
     143    iview_B=1;
     144end
     145RootPath_A=Param.InputTable{iview_A,1};
     146RootFile_A=Param.InputTable{iview_A,3};
     147SubDir_A=Param.InputTable{iview_A,2};
     148NomType_A=Param.InputTable{iview_A,4};
     149FileExt_A=Param.InputTable{iview_A,5};
     150RootPath_B=Param.InputTable{iview_B,1};
     151RootFile_B=Param.InputTable{iview_B,3};
     152SubDir_B=Param.InputTable{iview_B,2};
     153NomType_B=Param.InputTable{iview_B,4};
     154FileExt_B=Param.InputTable{iview_B,5};
     155
    154156
    155157PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
    156158
    157159[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,~,NomTypeNc]=...
    158                         find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
    159                    
     160    find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
     161
    160162if isempty(i1_series_Civ1)
    161163    disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun)
    … …  
    173175end
    174176Data.CivStage=0;%default
    175 list_param=(fieldnames(Param.ActionInput.Civ1))';
    176 list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
    177 Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
    178 Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
    179 Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
     177if Param.ActionInput.CheckCiv1
     178    list_param=(fieldnames(Param.ActionInput.Civ1))';
     179    list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
     180    Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
     181    Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
     182    Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
     183end
    180184% set the list of variables
    181185Data.ListVarName={'Coord_z','Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
    … …  
    224228    CheckOverwrite=Param.CheckOverwrite;
    225229end
    226    Data.Civ1_ImageA=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1,1));
    227    Data.Civ1_ImageB=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i1_series_Civ1(1),[],j2_series_Civ1(1,1));
    228     FileInfo=get_file_info(Data.Civ1_ImageA);
    229         par_civ1=Param.ActionInput.Civ1;% parameters for civ1
    230     par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height;
    231     par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width;
     230Data.Civ1_ImageA=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1,1));
     231Data.Civ1_ImageB=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i1_series_Civ1(1),[],j2_series_Civ1(1,1));
     232FileInfo=get_file_info(Data.Civ1_ImageA);
     233par_civ1=Param.ActionInput.Civ1;% parameters for civ1
     234par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height;
     235par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width;
    232236SearchRange_z=Param.ActionInput.Civ1.SearchRange(3);
    233     par_civ1.Dz=Param.ActionInput.Civ1.Dz;
    234     par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
     237par_civ1.Dz=Param.ActionInput.Civ1.Dz;
     238par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
    235239par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
    236240
    … …  
    268272        j2=j2_series_Civ1(ifield);
    269273    end
    270 
    271     Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair
    272     Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
    273 
    274     for ilist=1:length(list_param)
    275         Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
    276     end
    277    
    278     Data.CivStage=1;
    279  
     274   
     275   
     276   
    280277    ncfile_out=fullfile_uvmat(OutputPath,'',Param.InputTable{1,3},'.nc',...
    281                 '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
    282 
     278        '_1-1',i1_series_Civ1(ifield),i2_series_Civ1(ifield)); %output file name
     279   
    283280    if ~CheckOverwrite && exist(ncfile_out,'file')
    284281        disp(['existing output file ' ncfile_out ' already exists, skip to next field'])
    285282        continue% skip iteration if the mode overwrite is desactivated and the result file already exists
    286283    end
    287 
    288 
    289     %% Civ1
    290     % if Civ1 computation is requested
    291     if isfield (Param.ActionInput,'Civ1')
    292 
     284   
     285   
     286%% Civ1 computation if requested
     287    if Param.ActionInput.CheckCiv1
     288       
    293289        disp('civ1 started')
    294 
     290        Data.Civ1_Time=(Time(i2+1,j2+1)+Time(i1+1,j1+1))/2;% the Time is the Time at the middle of the image pair
     291        Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
     292       
     293        for ilist=1:length(list_param)
     294            Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
     295        end
     296       
     297        Data.CivStage=1;
    295298        % read input images by vertical blocks with nbre of images equal to 2*SearchRange+1,
    296299        par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);%image block initiation
    … …  
    307310            par_civ1.ImageB(iz,:,:) = B;
    308311        end
    309         % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
    310         [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
    311             Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), ~, errormsg] = civ3D (par_civ1);
    312         if ~isempty(errormsg)
    313             disp_uvmat('ERROR',errormsg,checkrun)
    314             return
    315         end
    316      % loop on slices
    317 
    318         for z_index=2:floor((NbSlice-SearchRange_z)/par_civ1.Dz)% loop on slices
    319             par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz,1);%shift the indices in the image block upward by par_civ1.Dz
    320             par_civ1.ImageB=circshift(par_civ1.ImageB,-par_civ1.Dz,1);
    321             for iz=1:par_civ1.Dz %read the new images at the end of the image block
    322                 image_index=z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz+1;
    323                 if image_index<=size(j1_series_Civ1,1)
    324                 j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz+1,1)
    325                 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%           
    326                 A= read_image(ImageName_A,FileType_A);
    327                 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
    328                 B= read_image(ImageName_B,FileType_B);
    329                 else
    330                     A=zeros(1,size(par_civ1.ImageA,2),size(par_civ1.ImageA,3));
    331                     B=zeros(1,size(par_civ1.ImageA,2),size(par_civ1.ImageA,3));
    332                 end
    333                 par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
    334                 par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
    335             end
    336             % caluclate velocity data (y and v in indices, reverse to y component)
    337             [Data.Civ1_X(z_index,:,:),Data.Civ1_Y(z_index,:,:), Data.Civ1_U(z_index,:,:), Data.Civ1_V(z_index,:,:),Data.Civ1_W(z_index,:,:),...
    338                 Data.Civ1_C(z_index,:,:), Data.Civ1_FF(z_index,:,:), result_conv, errormsg] = civ3D (par_civ1);
    339             if ~isempty(errormsg)
    340                 disp_uvmat('ERROR',errormsg,checkrun)
    341                 return
    342             end
    343         end
    344         Data.Civ1_V=-Data.Civ1_V;%reverse v
    345         Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
    346         [npz,npy,npx]=size(Data.Civ1_X);
    347        
    348       %  Data.Coord_y=flip(1:npy);
    349        % Data.Coord_x=1:npx;
     312       
    350313        % case of mask TO ADAPT
    351314        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
    … …  
    377340            end
    378341        end
     342       
     343        % calculate velocity data at the first position z, corresponding to j_index=par_civ1.Dz
     344        [Data.Civ1_X(1,:,:),Data.Civ1_Y(1,:,:), Data.Civ1_U(1,:,:), Data.Civ1_V(1,:,:),Data.Civ1_W(1,:,:),...
     345            Data.Civ1_C(1,:,:), Data.Civ1_FF(1,:,:), ~, errormsg] = civ3D (par_civ1);
     346        if ~isempty(errormsg)
     347            disp_uvmat('ERROR',errormsg,checkrun)
     348            return
     349        end
     350       
     351        % loop on slices
     352        for z_index=2:floor((NbSlice-SearchRange_z)/par_civ1.Dz)% loop on slices
     353            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz,1);%shift the indices in the image block upward by par_civ1.Dz
     354            par_civ1.ImageB=circshift(par_civ1.ImageB,-par_civ1.Dz,1);
     355            for iz=1:par_civ1.Dz %read the new images at the end of the image block
     356                image_index=z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz+1;
     357                if image_index<=size(j1_series_Civ1,1)
     358                    j_image_index=j1_series_Civ1(z_index*par_civ1.Dz+SearchRange_z-par_civ1.Dz+iz+1,1)
     359                    ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1,ifield),[],j_image_index);%
     360                    A= read_image(ImageName_A,FileType_A);
     361                    ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1,ifield),[],j_image_index);
     362                    B= read_image(ImageName_B,FileType_B);
     363                else
     364                    A=zeros(1,size(par_civ1.ImageA,2),size(par_civ1.ImageA,3));
     365                    B=zeros(1,size(par_civ1.ImageA,2),size(par_civ1.ImageA,3));
     366                end
     367                par_civ1.ImageA(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = A;
     368                par_civ1.ImageB(2*SearchRange_z+1-par_civ1.Dz+iz,:,:) = B;
     369            end
     370            % caluclate velocity data
     371            [Data.Civ1_X(z_index,:,:),Data.Civ1_Y(z_index,:,:), Data.Civ1_U(z_index,:,:), Data.Civ1_V(z_index,:,:),Data.Civ1_W(z_index,:,:),...
     372                Data.Civ1_C(z_index,:,:), Data.Civ1_FF(z_index,:,:), result_conv, errormsg] = civ3D (par_civ1);
     373            if ~isempty(errormsg)
     374                disp_uvmat('ERROR',errormsg,checkrun)
     375                return
     376            end
     377        end
     378        %         Data.Civ1_V=-Data.Civ1_V;%reverse v
     379        %         Data.Civ1_Y=npy-Data.Civ1_Y+1;%reverse y
     380       
     381       
     382        %  Data.Coord_y=flip(1:npy);
     383        % Data.Coord_x=1:npx;
     384       
     385    else
     386        Data=nc2struct(filecell{1,ifield});%read civ1 and fix1 data in the existing netcdf file
     387       
    379388        % Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
    380389        % Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
    381390        % Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];
    382         % 
    383         % 
     391        %
     392        %
    384393        % Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
    385394        % Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
    … …  
    389398        % Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
    390399        % Data.Civ1_FF=[Data.Civ1_FF reshape(F,[],1)];
    391 
    392     end
    393 
    394     %% Fix1
     400       
     401    end
     402   
     403%% Fix1
    395404    if isfield (Param.ActionInput,'Fix1')
    396405        disp('detect_false1 started')
    … …  
    412421        Data.CivStage=2;
    413422    end
    414     %% Patch1
     423   
     424%% Patch1
    415425    if isfield (Param.ActionInput,'Patch1')
    416426        disp('patch1 started')
    417427        tstart_patch1=tic;
    418 
     428       
    419429        % record the processing parameters of Patch1 as global attributes in the result nc file
    420430        list_param=fieldnames(Param.ActionInput.Patch1)';
    … …  
    426436        Data.CivStage=3;% record the new state of processing
    427437        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch1_param];
    428 
     438       
    429439        % list the variables to record
    430440        nbvar=length(Data.ListVarName);
    … …  
    433443        Data.VarAttribute{nbvar+1}.Role='vector_x';
    434444        Data.VarAttribute{nbvar+2}.Role='vector_y';
    435         Data.VarAttribute{nbvar+5}.Role='vector_z';     
     445        Data.VarAttribute{nbvar+5}.Role='vector_z';
    436446        Data.Civ1_U_smooth=Data.Civ1_U;
    437         Data.Civ1_V_smooth=Data.Civ1_V; 
    438         Data.Civ1_W_smooth=Data.Civ1_W; 
     447        Data.Civ1_V_smooth=Data.Civ1_V;
     448        Data.Civ1_W_smooth=Data.Civ1_W;
    439449        if isfield(Data,'Civ1_FF')
    440450            ind_good=find(Data.Civ1_FF==0);
    … …  
    447457            return
    448458        end
    449 
     459       
    450460        % perform Patch calculation using the UVMAT fct 'filter_tps'
    451 Civ1_Z=0.5*Data.Civ1_W(ind_good);
    452 for iz=1:numel(Data.Coord_z)
    453     Civ1_Z(iz,:,:)=Civ1_Z(iz,:,:)+Data.Coord_z(iz);
    454 end
     461        Civ1_Z=0.5*Data.Civ1_W;
     462        for iz=1:numel(Data.Coord_z)
     463            Civ1_Z(iz,:,:)=Civ1_Z(iz,:,:)+Data.Coord_z(iz);   
     464        end
     465        [npz,npy,npx]=size(Data.Civ1_X);
    455466        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,Data.Civ1_W_tps,...
    456             Data.Civ1_U_smooth(ind_good),Data.Civ1_V_smooth(ind_good),Data.Civ1_V_smooth(ind_good),FFres]=...
    457             filter_tps_3D([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Civ_1_Z,Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),Data.Civ1_W(ind_good),Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
    458         Data.Civ1_FF(ind_good)=uint8(FFres);
     467            Data.Civ1_U_smooth(ind_good),Data.Civ1_V_smooth(ind_good),Data.Civ1_W_smooth(ind_good),FFres]=...
     468            filter_tps_3D(Data.Civ1_X(ind_good),Data.Civ1_Y(ind_good),Civ1_Z(ind_good),Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),Data.Civ1_W(ind_good),...
     469            Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
     470        Data.Civ1_FF(ind_good)=uint8(4*FFres);
     471        Data.Civ1_FF=reshape(Data.Civ1_FF,npz,npy,npx);
     472        Data.Civ1_U=reshape(Data.Civ1_U,npz,npy,npx);
     473        Data.Civ1_V=reshape(Data.Civ1_V,npz,npy,npx);
     474        Data.Civ1_W=reshape(Data.Civ1_W,npz,npy,npx);
    459475        time_patch1=toc(tstart_patch1);
    460476        disp('patch1 performed')
    461477    end
    462 
     478   
    463479    %% Civ2
    464480    if isfield (Param.ActionInput,'Civ2')
    … …  
    503519        end
    504520        % end
    505 
     521       
    506522        % get the guess from patch1 or patch2 (case 'CheckCiv3')
    507 
     523       
    508524        if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
    509525            SubRange= Data.Civ2_SubRange;
    … …  
    530546            Data.CivStage=4;
    531547        end
    532 
     548       
    533549        Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
    534550        Shifty=zeros(size(par_civ2.Grid,1),1);
    535551        nbval=zeros(size(par_civ2.Grid,1),1);% nbre of interpolated values at each grid point (from the different patch subdomains)
    536 
     552       
    537553        NbSubDomain=size(SubRange,3);
    538554        for isub=1:NbSubDomain% for each sub-domain of Patch1
    … …  
    577593            end
    578594        end
    579 
    580 
     595       
     596       
    581597        if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    582598            Civ1_Dt=1;
    … …  
    585601            Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1);
    586602        end
    587 
     603       
    588604        par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
    589605        % shift the grid points by half the expected shift to provide the correlation box position in image A
    … …  
    595611            par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
    596612        end
    597 
     613       
    598614        % calculate velocity data (y and v in image indices, reverse to y component)
    599 
     615       
    600616        [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2);
    601 
     617       
    602618        list_param=(fieldnames(Param.ActionInput.Civ2))';
    603619        list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list
    … …  
    620636        end
    621637        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param];
    622 
     638       
    623639        nbvar=numel(Data.ListVarName);
    624640        % define the Civ2 variable (if Civ2 data are not replaced from previous calculation)
    … …  
    651667        end
    652668    end
    653 
     669   
    654670    %% Fix2
    655671    if isfield (Param.ActionInput,'Fix2')
    … …  
    665681        Data.CivStage=Data.CivStage+1;
    666682    end
    667 
     683   
    668684    %% Patch2
    669685    if isfield (Param.ActionInput,'Patch2')
    … …  
    678694        end
    679695        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch2_param];
    680 
     696       
    681697        nbvar=length(Data.ListVarName);
    682698        Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
    683699        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
    684700            {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
    685 
     701       
    686702        Data.VarAttribute{nbvar+1}.Role='vector_x';
    687703        Data.VarAttribute{nbvar+2}.Role='vector_y';
    … …  
    700716            return
    701717        end
    702 
     718       
    703719        [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=...
    704720            filter_tps([Data.Civ2_X(ind_good) Data.Civ2_Y(ind_good)],Data.Civ2_U(ind_good),Data.Civ2_V(ind_good),[],Data.Patch2_SubDomainSize,Data.Patch2_FieldSmooth,Data.Patch2_MaxDiff);
    … …  
    710726        disp('patch2 performed')
    711727    end
    712 
     728   
    713729    %% write result in a netcdf file if requested
    714730    % if CheckOutputFile
    … …  
    726742    disp(['time civ2 ' num2str(time_civ2,2) ' s'])
    727743    disp(['time patch2 ' num2str(time_patch2,2) ' s'])
    728 
     744   
    729745end
    730746
    … …  
    747763%  .ImageB: second image for correlation(matrix)
    748764%  .CorrBoxSize: 1,2 vector giving the size of the correlation box in x and y
    749 %  .SearchBoxSize:  1,2 vector giving the size of the search box in x and y
     765%  .SearchRange:  1,2 vector giving the range of the search in x and y
    750766%  .SearchBoxShift: 1,2 vector or 2 column matrix (for civ2) giving the shift of the search box in x and y
    751767%  .CorrSmooth: =1 or 2 determines the choice of the sub-pixel determination of the correlation max
    … …  
    764780function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ)
    765781%% check image sizes
    766 [npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz = 
     782[npz,npy_ima,npx_ima]=size(par_civ.ImageA);% npz =
    767783if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
    768784    errormsg='image pair with unequal size';
    … …  
    786802ibx2=floor(par_civ.CorrBoxSize(1)/2);
    787803iby2=floor(par_civ.CorrBoxSize(2)/2);
    788 isx2=par_civ.SearchRange(1);
    789 isy2=par_civ.SearchRange(2);
     804isx2=par_civ.SearchRange(1)+ibx2;
     805isy2=par_civ.SearchRange(2)+iby2;
    790806isz2=par_civ.SearchRange(3);
    791807kref=isz2+1;%middle index of the z slice
    … …  
    878894                image2_mean=mean(mean(image2_crop,2),3);
    879895            end
    880      
     896            
    881897            if FF(ivec_y,ivec_x)==0
    882898                if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    883899                    FF(ivec_y,ivec_x)=1;       %threshold on image minimum
    884                
     900                   
    885901                elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
    886902                    FF(ivec_y,ivec_x)=1;    %threshold on image maximum
    887903                end
    888904                if FF(ivec_y,ivec_x)==1
    889                     utable(ivec_y,ivec_x)=NaN; 
     905                    utable(ivec_y,ivec_x)=NaN;
    890906                    vtable(ivec_y,ivec_x)=NaN;
    891907                else
    … …  
    898914                    %     % image2_crop=(image2_crop-image2_mean);
    899915                    % end
    900 
     916                   
    901917                    npxycorr=2*par_civ.SearchRange(1:2)+1;
    902918                    result_conv=zeros([2*par_civ.SearchRange(3)+1 npxycorr]);%initialise the conv product
    903919                    max_xy=zeros(2*par_civ.SearchRange(3)+1,1);%initialise the max correlation vs z
    904                     xk=ones(npz,1);%initialise the x displacement corresponding to the max corr 
    905                     yk=ones(npz,1);%initialise the y displacement corresponding to the max corr 
     920                    xk=ones(npz,1);%initialise the x displacement corresponding to the max corr
     921                    yk=ones(npz,1);%initialise the y displacement corresponding to the max corr
    906922                    subima1=squeeze(image1_crop(kref,:,:))-image1_mean(kref);
    907923                    for kz=1:npz
    908                         subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz);   
     924                        subima2=squeeze(image2_crop(kz,:,:))-image2_mean(kz);
    909925                        correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid');
    910926                        result_conv(kz,:,:)=correl_xy;
    911                         max_xy(kz)=max(max(correl_xy));             
    912                          [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);%find the indices of the maximum, use 0.999 to avoid rounding errors
     927                        max_xy(kz)=max(max(correl_xy));
     928                        [yk(kz),xk(kz)]=find(correl_xy==max_xy(kz),1);%find the indices of the maximum, use 0.999 to avoid rounding errors
    913929                    end
    914930                    [corrmax,dz]=max(max_xy);
    … …  
    932948                        end
    933949                        utable(ivec_y,ivec_x)=vector(1)+shiftx(ivec_y,ivec_x);
    934                         vtable(ivec_y,ivec_x)=vector(2)+shifty(ivec_y,ivec_x);
     950                        vtable(ivec_y,ivec_x)=-vector(2)+shifty(ivec_y,ivec_x);
    935951                        wtable(ivec_y,ivec_x)=vector(3);
    936952                        xtable(ivec_y,ivec_x)=iref+utable(ivec_y,ivec_x)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    … …  
    938954                        iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
    939955                        jref=round(ytable(ivec_y,ivec_x)+0.5);
    940                        % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0
    941                        %     'test'
    942                        % end
     956                        % if utable(ivec_y,ivec_x)==0 && vtable(ivec_y,ivec_x)==0
     957                        %     'test'
     958                        % end
    943959                        % eliminate vectors located in the mask
    944960                        if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    … …  
    956972    end
    957973end
     974ytable=npy_ima-ytable+1;%reverse from j index to image coordinate y
    958975result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
    959976
    … …  
    980997    f2 = log(result_conv(z+1,y,x));
    981998    peakz = peakz+ (f1-f2)/(2*f1-4*f0+2*f2);
    982 
     999   
    9831000    f1 = log(result_conv(z,y-1,x));
    9841001    f2 = log(result_conv(z,y+1,x));
    9851002    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    986 
     1003   
    9871004    f1 = log(result_conv(z,y,x-1));
    9881005    f2 = log(result_conv(z,y,x+1));
    … …  
    10841101
    10851102if isfield (Param,'MinCorr')
    1086      FF(C<Param.MinCorr & FFIn==0)=2;
     1103    FF(C<Param.MinCorr & FFIn==0)=2;
    10871104end
    10881105if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
    … …  
    10931110    end
    10941111    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
    1095          U2Max=Param.MaxVel*Param.MaxVel;
     1112        U2Max=Param.MaxVel*Param.MaxVel;
    10961113        FF(Umod>U2Max & FFIn==0)=3;
    10971114    end
    … …  
    11391156        i1_series=i_series-ind1;% set of first image numbers
    11401157        i2_series=i_series+ind2;
    1141         check_bounds=i1_series<MinIndex_i | i2_series>MaxIndex_i;
     1158        check_bounds=i1_series(1)<MinIndex_i | i2_series(1)>MaxIndex_i;
    11421159        if isempty(j_series)
    11431160            NomTypeNc='_1-2';
Note: See TracChangeset for help on using the changeset viewer.