Ignore:
Timestamp:
Jun 21, 2024, 4:51:59 PM (3 months ago)
Author:
sommeria
Message:

civ_3D corrected

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/civ_3D.m

    r1127 r1148  
    1 %'civ_3D': PIV function activated by the general GUI series
    2 % --- call the sub-functions:
    3 %   civ: PIV function itself
    4 %   fix: removes false vectors after detection by various criteria
    5 %   filter_tps: make interpolation-smoothing
     1%'civ_3D': 3D PIV from image scan in a volume
     2
    63%------------------------------------------------------------------------
    7 % function [Data,errormsg,result_conv]= civ_series(Param)
     4% function [Data,errormsg,result_conv]= civ_3D(Param)
    85%
    96%OUTPUT
    … …  
    2118%           if absent no file is produced, result in the output structure Data (test mode)
    2219%     Param.ActionInput: substructure with the parameters provided by the GUI civ_input
    23 %                      .Civ1: parameters for civ1
    24 %                      .Fix1: parameters for fix1
     20%                      .Civ1: parameters for civ1cc
     21%                      .Fix1: parameters for detect_false1
    2522%                      .Patch1:
    2623%                      .Civ2: for civ2
    … …  
    5350    path_series=fileparts(which('series'));
    5451    addpath(fullfile(path_series,'series'))
    55     Data=civ_input_3D(Param);% introduce the civ parameters using the GUI civ_input
    56     %Data=civ_input_App(Param);% introduce the civ parameters using the GUI
    57     %civ_input %%%%  TO MOVE
    58     if isempty(Data)
    59         Data=Param;% if  civ_input has been cancelled, keep previous parameters
    60     end
     52
     53   Data=civ_input(Param)
     54
    6155    Data.Program=mfilename;%gives the name of the current function
    6256    Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    63     Data.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
     57    Data.WholeIndexRange_j='on';% prescribes the file index ranges j from min to max (options 'off'/'on', 'off' by default)
    6458    Data.NbSlice='off'; %nbre of slices ('off' by default)
    6559    Data.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    … …  
    6862    Data.ProjObject='off';%can use projection object(option 'off'/'on',
    6963    Data.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
    70     Data.OutputDirExt='.civ';%set the output dir extension
     64    Data.OutputDirExt='.civ_3D';%set the output dir extension
    7165    Data.OutputSubDirMode='last'; %select the last subDir in the input table as root of the output subdir name (option 'all'/'first'/'last', 'all' by default)
    7266    Data.OutputFileMode='NbInput_i';% one output file expected per value of i index (used for waitbar)
    7367    Data.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1)
    74     if isfield(Data.ActionInput,'PairIndices') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
     68    if isfield(Data,'ActionInput') && isfield(Data.ActionInput,'PairIndices') && strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
    7569        if isfield(Data.ActionInput.PairIndices,'ListPairCiv2')
    7670            str_civ=Data.ActionInput.PairIndices.ListPairCiv2;
    … …  
    9084        end
    9185    end
    92     % estimated CPUTime
    93     CPUtime_unit=0.01;%estimated time for a multiplication (in microsecond)
    94     if isfield(Param.SeriesData,'FileInfo')&&isfield(Param.SeriesData.FileInfo{1},'Height')&&isfield(Param.SeriesData.FileInfo{1},'Width')
    95         pixnbre=Param.SeriesData.FileInfo{1}.Height*Param.SeriesData.FileInfo{1}.Width; % total number of pxels for input images 
    96         CPUtime=0;
    97         if isfield(Data.ActionInput,'Civ1')
    98             %BoxSize=Data.ActionInput.Civ1.CorrBoxSize(1)*Data.ActionInput.Civ1.CorrBoxSize(2);
    99             tic
    100             testboxa=rand(Data.ActionInput.Civ1.CorrBoxSize(1),Data.ActionInput.Civ1.CorrBoxSize(2));
    101             testboxb=rand(Data.ActionInput.Civ1.SearchBoxSize(1),Data.ActionInput.Civ1.SearchBoxSize(2));
    102             anss=conv2(testboxa,testboxb);
    103             CPUtime_unit=toc;
    104             nb_box=pixnbre/(Data.ActionInput.Civ1.Dx*Data.ActionInput.Civ1.Dy);   
    105             %nbpos=Data.ActionInput.Civ1.SearchBoxSize-Data.ActionInput.Civ1.CorrBoxSize;
    106             CPUtime=2*CPUtime_unit*nb_box%*BoxSize*nbpos(1)*nbpos(2);% adjustement factor 2 used
    107         end
    108         if isfield(Data.ActionInput,'Patch1')
    109             CPUtime=2*CPUtime;
    110         end
    111         if isfield(Data.ActionInput,'Civ2')
    112             tic
    113             testboxa=rand(Data.ActionInput.Civ2.CorrBoxSize(1),Data.ActionInput.Civ2.CorrBoxSize(2));
    114             testboxb=rand(Data.ActionInput.Civ2.SearchBoxSize(1),Data.ActionInput.Civ2.SearchBoxSize(2));
    115             anss=conv2(testboxa,testboxb);
    116             CPUtime_unit=toc;
    117             nb_box=pixnbre/(Data.ActionInput.Civ2.Dx*Data.ActionInput.Civ2.Dy);
    118             %BoxSize=Data.ActionInput.Civ2.CorrBoxSize(1)*Data.ActionInput.Civ2.CorrBoxSize(2);
    119             %nbpos=Data.ActionInput.Civ2.SearchBoxSize-Data.ActionInput.Civ2.CorrBoxSize;
    120             CPUtime=CPUtime+2*CPUtime_unit*nb_box;%*BoxSize*nbpos(1)*nbpos(2);
    121         end
    122         if isfield(Data.ActionInput,'Patch2')
    123             CPUtime=(4/3)*CPUtime;
    124         end
    125         Data.CPUTime=ceil(CPUtime/6); % estimated CPU time per field pair in minute
    126         Data.CPUTime=Data.CPUTime/10; % displqy CPU time with 1 digit beyond dot
    127     end
     86
     87
    12888    return
    12989end
    130 
    131 %% read input parameters from an xml file if input is a file name (batch mode)
    132 checkrun=1;
     90%% END OF ENTERING INPUT PARAMETER MODE
     91
     92%% RUN MODE: read input parameters from an xml file if input is a file name (batch mode)
    13393if ischar(Param)
    13494    Param=xml2struct(Param);% read Param as input file (batch case)
    13595    checkrun=0;
    136 end
    137 
    138 %% test input
     96    RUNHandle=[];
     97else
     98    hseries=findobj(allchild(0),'Tag','series');
     99    RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
     100    WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
     101    checkrun=1;
     102end
     103
     104%% test input Param
     105if ~isfield(Param,'InputTable')
     106    disp('ERROR: no input file entered')
     107    return
     108end
    139109if ~isfield(Param,'ActionInput')
    140110    disp_uvmat('ERROR','no parameter set for PIV',checkrun)
    … …  
    142112end
    143113iview_A=0;%default values
    144 NbField=1;
    145 RUNHandle=[];
    146 CheckInputFile=isfield(Param,'InputTable');%= 1 in test use for TestCiv (no nc file involved)
    147 CheckOutputFile=isfield(Param,'OutputSubDir');%= 1 in test use for TestPatch (no nc file produced)
    148 
    149 %% input files and indexing (skipped in Test mode)
    150 if CheckInputFile
    151     hseries=findobj(allchild(0),'Tag','series');
    152     RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
    153     WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
    154     MaxIndex_i=Param.IndexRange.MaxIndex_i;
    155     MinIndex_i=Param.IndexRange.MinIndex_i;
    156     MaxIndex_j=ones(size(MaxIndex_i));MinIndex_j=ones(size(MinIndex_i));
    157     if isfield(Param.IndexRange,'MaxIndex_j')&& isfield(Param.IndexRange,'MinIndex_j')
    158         MaxIndex_j=Param.IndexRange.MaxIndex_j;
    159         MinIndex_j=Param.IndexRange.MinIndex_j;
    160     end
    161     if isfield(Param,'InputTable')
    162         [tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
    163         iview_A=0;% series index (iview) for the first image series
    164         iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
    165         if Param.ActionInput.CheckCiv1
    166             iview_A=1;% usual PIV, the image series is on the first line of the table
    167         elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line
    168             iview_A=2;% the second line is used for the input images of Civ2
    169         end
    170 %         if strcmp(Param.ActionInput.ListCompareMode,'shift')
    171 %             iview_B=iview_A+1; % the second image series is on the next line of the input table
    172 %         end
    173         if iview_A~=0
    174             RootPath_A=Param.InputTable{iview_A,1};
    175             RootFile_A=Param.InputTable{iview_A,3};
    176             SubDir_A=Param.InputTable{iview_A,2};
    177             NomType_A=Param.InputTable{iview_A,4};
    178             FileExt_A=Param.InputTable{iview_A,5};
    179             if iview_B==0
    180                 iview_B=iview_A;% the second image series is the same as the first
    181             end
    182             RootPath_B=Param.InputTable{iview_B,1};
    183             RootFile_B=Param.InputTable{iview_B,3};
    184             SubDir_B=Param.InputTable{iview_B,2};
    185             NomType_B=Param.InputTable{iview_B,4};
    186             FileExt_B=Param.InputTable{iview_B,5};
    187         end
    188        
    189         PairCiv2='';
    190         switch Param.ActionInput.ListCompareMode
    191             case 'PIV'
    192                 PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
    193                 if isfield(Param.ActionInput.PairIndices,'ListPairCiv2')
    194                     PairCiv2=Param.ActionInput.PairIndices.ListPairCiv2;%string which determines the civ2 pair
    195                 end
    196                 if iview_A==1% if Civ1 is performed
    197                     [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
     114
     115if isfield(Param,'OutputSubDir')&& isfield(Param,'OutputDirExt')
     116    OutputDir=[Param.OutputSubDir Param.OutputDirExt];
     117     OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
     118else
     119    disp_uvmat('ERROR','no output folder defined',checkrun)
     120    return
     121end
     122
     123%% input files and indexing
     124MaxIndex_i=Param.IndexRange.MaxIndex_i;
     125MinIndex_i=Param.IndexRange.MinIndex_i;
     126MaxIndex_j=ones(size(MaxIndex_i));MinIndex_j=ones(size(MinIndex_i));
     127if isfield(Param.IndexRange,'MaxIndex_j')&& isfield(Param.IndexRange,'MinIndex_j')
     128    MaxIndex_j=Param.IndexRange.MaxIndex_j;
     129    MinIndex_j=Param.IndexRange.MinIndex_j;
     130end
     131
     132[tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     133iview_A=0;% series index (iview) for the first image series
     134iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
     135if Param.ActionInput.CheckCiv1
     136    iview_A=1;% usual PIV, the image series is on the first line of the table
     137elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line
     138    iview_A=2;% the second line is used for the input images of Civ2
     139end
     140if iview_A~=0
     141    RootPath_A=Param.InputTable{iview_A,1};
     142    RootFile_A=Param.InputTable{iview_A,3};
     143    SubDir_A=Param.InputTable{iview_A,2};
     144    NomType_A=Param.InputTable{iview_A,4};
     145    FileExt_A=Param.InputTable{iview_A,5};
     146    if iview_B==0
     147        iview_B=iview_A;% the second image series is the same as the first
     148    end
     149    RootPath_B=Param.InputTable{iview_B,1};
     150    RootFile_B=Param.InputTable{iview_B,3};
     151    SubDir_B=Param.InputTable{iview_B,2};
     152    NomType_B=Param.InputTable{iview_B,4};
     153    FileExt_B=Param.InputTable{iview_B,5};
     154end
     155
     156PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
     157
     158[i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
    198159                        find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
    199                     if ~isempty(PairCiv2)
    200                         [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds_Civ2]=...
    201                             find_pair_indices(PairCiv2,i1_series{1},j1_series{1},MinIndex_i(1),MaxIndex_i(1),MinIndex_j(1),MaxIndex_j(1));
    202                         check_bounds=check_bounds | check_bounds_Civ2;
    203                     end
    204                 else% we start from an existing Civ1 file
    205                     i1_series_Civ1=i1_series{1};
    206                     i2_series_Civ1=i2_series{1};
    207                     j1_series_Civ1=j1_series{1};
    208                     j2_series_Civ1=j2_series{1};
    209                     NomTypeNc=Param.InputTable{1,4};
    210                     if ~isempty(PairCiv2)
    211                         [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds,NomTypeNc]=...
    212                             find_pair_indices(PairCiv2,i1_series{2},j1_series{2},MinIndex_i(2),MaxIndex_i(2),MinIndex_j(2),MaxIndex_j(2));
    213                     end
    214                 end
    215             case 'displacement'
    216                 if isfield(Param.ActionInput,'OriginIndex')
    217                 i1_series_Civ1=Param.ActionInput.OriginIndex*ones(size(i1_series{1}));
    218                 else
    219                     i1_series_Civ1=ones(size(i1_series{1}));
    220                 end
    221                 i1_series_Civ2=i1_series_Civ1;
    222                 i2_series_Civ1=i1_series{1};
    223                 i2_series_Civ2=i1_series{1};
    224                 j1_series_Civ1=[];% no j index variation for the ref image
    225                 j1_series_Civ2=[];
    226                 if isempty(j1_series{1})
    227                     j2_series_Civ1=ones(size(i1_series_Civ1));
    228                 else
    229                     j2_series_Civ1=j1_series{1};% if j index exist 
    230                 end
    231                 j2_series_Civ2=j2_series_Civ1;
    232                 NomTypeNc='_1';
    233             case 'PIV volume'
    234                 % TODO, TODO
    235         end
    236         %determine frame indices for input with movie or other multiframe input file
    237         if isempty(j1_series_Civ1)% simple movie with index i
    238             FrameIndex_A_Civ1=i1_series_Civ1;
    239             FrameIndex_B_Civ1=i2_series_Civ1;
    240             j1_series_Civ1=ones(size(i1_series_Civ1));
    241             if strcmp(Param.ActionInput.ListCompareMode,'PIV')
    242             j2_series_Civ1=ones(size(i1_series_Civ1));
    243             end
    244         else % movie for each burst or volume (index j)
    245             FrameIndex_A_Civ1=j1_series_Civ1;
    246             FrameIndex_B_Civ1=j2_series_Civ1;
    247         end
    248         if isempty(PairCiv2)
    249             FrameIndex_A_Civ2=FrameIndex_A_Civ1;
    250             FrameIndex_B_Civ2=FrameIndex_B_Civ1;
    251         else
    252             if isempty(j1_series_Civ2)
    253                 FrameIndex_A_Civ2=i1_series_Civ2;
    254                 FrameIndex_B_Civ2=i2_series_Civ2;
    255                 j1_series_Civ2=ones(size(i1_series_Civ2));
    256                 if strcmp(Param.ActionInput.ListCompareMode,'PIV')
    257                 j2_series_Civ2=ones(size(i1_series_Civ2));
    258                 end
    259             else
    260                 FrameIndex_A_Civ2=j1_series_Civ2;
    261                 FrameIndex_B_Civ2=j2_series_Civ2;
    262             end
    263         end
    264         if isempty(i1_series_Civ1)||(~isempty(PairCiv2) && isempty(i1_series_Civ2))
    265             disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun)
    266             return
    267         end
    268     end
    269    
    270     %% check the first image pair
    271         if Param.ActionInput.CheckCiv1% Civ1 is performed
    272             NbField=numel(i1_series_Civ1);
    273         elseif Param.ActionInput.CheckCiv2 % Civ2 is performed without Civ1
    274             NbField=numel(i1_series_Civ2);
    275         else
    276             NbField=numel(i1_series_Civ1);% no image used (only fix or patch) TO CHECK
    277         end
    278 
    279     %% Output directory
    280     OutputDir='';
    281     if CheckOutputFile
    282         OutputDir=[Param.OutputSubDir Param.OutputDirExt];
    283     end
    284 end
     160                   
     161if isempty(i1_series_Civ1)
     162    disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun)
     163    return
     164end
     165NbField_i=size(i1_series_Civ1,2);
     166NbSlice=size(i1_series_Civ1,1);
    285167
    286168%% prepare output Data
    … …  
    289171Data.Program='civ_series';
    290172Data.CivStage=0;%default
    291 check_civx=0;%default
     173list_param=(fieldnames(Param.ActionInput.Civ1))';
     174    list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
     175    Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
     176    Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
     177Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
     178            % set the list of variables
     179        Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_FF'};%  cell array containing the names of the fields to record
     180        Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
     181        Data.VarAttribute{1}.Role='coord_x';
     182        Data.VarAttribute{2}.Role='coord_y';
     183        Data.VarAttribute{3}.Role='vector_x';
     184        Data.VarAttribute{4}.Role='vector_y';
     185        Data.VarAttribute{5}.Role='vector_z';
     186        Data.VarAttribute{6}.Role='ancillary';
     187        Data.VarAttribute{7}.Role='errorflag';
    292188
    293189%% get timing from the ImaDoc file or input video
    … …  
    314210end
    315211
    316 %%%%% MAIN LOOP %%%%%%
    317212maskoldname='';% initiate the mask name
    318213FileType_A='';
    … …  
    322217    CheckOverwrite=Param.CheckOverwrite;
    323218end
    324 for ifield=1:NbField
    325     tic
    326     if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
    327         update_waitbar(WaitbarHandle,ifield/NbField)
     219   Data.Civ1_ImageA=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1,1));
     220   Data.Civ1_ImageB=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i1_series_Civ1(1),[],j2_series_Civ1(1,1));
     221    FileInfo=get_file_info(Data.Civ1_ImageA);
     222        par_civ1=Param.ActionInput.Civ1;% parameters for civ1
     223    par_civ1.ImageHeight=FileInfo.Height;npy=FileInfo.Height;
     224    par_civ1.ImageWidth=FileInfo.Width;npx=FileInfo.Width;
     225SearchRange_z=floor(Param.ActionInput.Civ1.SearchBoxSize(3)/2);
     226    par_civ1.Dz=Param.ActionInput.Civ1.Dz;
     227    par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
     228par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
     229
     230%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     231%%%%% MAIN LOOP %%%%%%
     232%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     233for ifield=1:NbField_i
     234    tstart=tic;
     235    time_civ1=0;
     236    time_patch1=0;
     237    time_civ2=0;
     238    time_patch2=0;
     239    if checkrun% update the waitbar in interactive mode with GUI series  (checkrun=1)
     240        update_waitbar(WaitbarHandle,ifield/NbField_i)
    328241        if  checkrun && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    329242            disp('program stopped by user')
    … …  
    331244        end
    332245    end
    333     if CheckInputFile
    334         if iview_A==0 % no nc file has been entered
    335             ncfile=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...
    336                 NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));
    337         else% an existing nc file has been entered
    338             if iview_A==1% if Civ1 is performed
    339                 Civ1Dir=OutputDir;
    340             else
    341                 Civ1Dir=Param.InputTable{1,2};
    342             end
    343             if strcmp(Param.ActionInput.ListCompareMode,'PIV')
    344                 ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),...
    345                     j1_series_Civ1(ifield),j2_series_Civ1(ifield));
    346             else
    347                 ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
    348                     j1_series_Civ1(ifield),j2_series_Civ1(ifield));
    349             end
    350         end
    351         ncfile_out=ncfile;% by default
    352         if isfield (Param.ActionInput,'Civ2')
    353             i1_civ2=i1_series_Civ2(ifield);
    354             i2_civ2=i1_civ2;
    355             if ~isempty(i2_series_Civ2)
    356                 i2_civ2=i2_series_Civ2(ifield);
    357             end
    358             j1_civ2=1;
    359             if ~isempty(j1_series_Civ2)
    360                 j1_civ2=j1_series_Civ2(ifield);
    361             end
    362             j2_civ2=i1_civ2;
    363             if ~isempty(j2_series_Civ2)
    364                 j2_civ2=j2_series_Civ2(ifield);
    365             end
    366             if strcmp(Param.ActionInput.ListCompareMode,'PIV')
    367                 ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1_civ2,i2_civ2,j1_civ2,j2_civ2);
    368             else % displacement
    369                 ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i2_civ2,[],j2_civ2);
    370             end
    371         end
    372         if ~CheckOverwrite && exist(ncfile_out,'file')
    373             disp(['existing output file ' ncfile_out ' already exists, skip to next field'])
    374             continue% skip iteration if the mode overwrite is desactivated and the result file already exists
    375         end
    376     end
     246   
     247    %indicate the values of all the global attributes in the output data
     248    i1=i1_series_Civ1(ifield);
     249    i2=i1;
     250    if ~isempty(i2_series_Civ1)
     251        i2=i2_series_Civ1(ifield);
     252    end
     253    j1=1;
     254    if ~isempty(j1_series_Civ1)
     255        j1=j1_series_Civ1(ifield);
     256    end
     257    j2=j1;
     258    if ~isempty(j2_series_Civ1)
     259        j2=j2_series_Civ1(ifield);
     260    end
     261
     262    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
     263    Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
     264
     265    for ilist=1:length(list_param)
     266        Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
     267    end
     268   
     269    Data.CivStage=1;
     270    % output nc file
     271    ncfile_out=fullfile_uvmat(OutputPath,Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...
     272        NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));
     273
     274    if ~CheckOverwrite && exist(ncfile_out,'file')
     275        disp(['existing output file ' ncfile_out ' already exists, skip to next field'])
     276        continue% skip iteration if the mode overwrite is desactivated and the result file already exists
     277    end
     278
     279
    377280    %% Civ1
    378281    % if Civ1 computation is requested
    379282    if isfield (Param.ActionInput,'Civ1')
    380         if CheckInputFile
    381             disp('civ1 started')
    382         end
    383         par_civ1=Param.ActionInput.Civ1;% parameters for civ1
    384         if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
    385             try
    386                 if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    387                     ImageName_A=Param.ActionInput.RefFile;
    388                 else
    389                 ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
    390                 end
    391                 if strcmp(FileExt_A,'.nc')% case of input images in format netcdf
    392                     FieldName_A=Param.InputFields.FieldName;
    393                     [DataIn,tild,tild,errormsg]=nc2struct(ImageName_A,{FieldName_A});
    394                     par_civ1.ImageA=DataIn.(FieldName_A);
    395                 else % usual image formats for image A
    396                     if isempty(FileType_A)% open the image object if not already done in case of movie input
    397                         [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A);
    398                         FileType_A=FileInfo_A.FileType;
    399                         if isempty(Time) && ~isempty(find(strcmp(FileType_A,{'mmreader','video','cine_phantom'})))% case of video input
    400                             Time=zeros(FileInfo_A.NumberOfFrames+1,2);
    401                             Time(:,2)=(0:1/FileInfo_A.FrameRate:(FileInfo_A.NumberOfFrames)/FileInfo_A.FrameRate)';
    402                         end
    403                         if ~isempty(FileType_A) && isempty(Time)% Time = index i +0.001 index j by default
    404                             MaxIndex_i=max(i2_series_Civ1);
    405                             MaxIndex_j=max(j2_series_Civ1);
    406                             Time=(1:MaxIndex_i)'*ones(1,MaxIndex_j);
    407                             Time=Time+0.001*ones(MaxIndex_i,1)*(1:MaxIndex_j);
    408                             Time=[zeros(1,MaxIndex_j);Time];% insert a first line of zeros
    409                             Time=[zeros(MaxIndex_i+1,1) Time];% insert a first column of zeros
    410                         end
    411                     end
    412                     if ~exist(ImageName_A,'file')
    413                         disp([ImageName_A ' missing'])
    414                         continue
    415                     end
    416                     [par_civ1.ImageA,VideoObject_A] = read_image(ImageName_A,FileType_A,VideoObject_A,FrameIndex_A_Civ1(ifield));
    417                 end
    418                 ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield));
    419                 if strcmp(FileExt_B,'.nc') % case of input images in format netcdf
    420                     FieldName_B=Param.InputFields.FieldName;
    421                     [DataIn,tild,tild,errormsg]=nc2struct(ImageName_B,{FieldName_B});
    422                     par_civ1.ImageB=DataIn.(FieldName_B);
    423                 else % usual image formats for image B
    424                     if isempty(FileType_B)
    425                         [FileInfo_B,VideoObject_B]=get_file_info(ImageName_B);
    426                         FileType_B=FileInfo_B.FileType;
    427                     end
    428                     if ~exist(ImageName_B,'file')
    429                         disp([ImageName_B ' missing'])
    430                         continue
    431                     end
    432                     [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
    433                 end
    434             catch ME % display errors in reading input images
    435                 if ~isempty(ME.message)
    436                     disp_uvmat('ERROR', ['error reading input image: ' ME.message],checkrun)
    437                     continue
    438                 end
    439             end
    440             par_civ1.ImageWidth=size(par_civ1.ImageA,2);
    441             par_civ1.ImageHeight=size(par_civ1.ImageA,1);
    442             list_param=(fieldnames(Param.ActionInput.Civ1))';
    443             list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
    444             Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
    445             Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
    446             %indicate the values of all the global attributes in the output data
    447             Data.Civ1_ImageA=ImageName_A;
    448             Data.Civ1_ImageB=ImageName_B;
    449             i1=i1_series_Civ1(ifield);
    450             i2=i1;
    451             if ~isempty(i2_series_Civ1)
    452                 i2=i2_series_Civ1(ifield);
    453             end
    454             j1=1;
    455             if ~isempty(j1_series_Civ1)
    456                 j1=j1_series_Civ1(ifield);
    457             end
    458             j2=j1;
    459             if ~isempty(j2_series_Civ1)
    460                 j2=j2_series_Civ1(ifield);
    461             end
    462             if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    463                 Data.Civ1_Time=Time(i2+1,j2+1);% the Time is the Time of the secodn image
    464                 Data.Civ1_Dt=1;% Time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
    465             else
    466                 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
    467                 Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
    468             end
    469             for ilist=1:length(list_param)
    470                 Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
    471             end
    472             Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
    473             Data.CivStage=1;
    474         end
    475         % set the list of variables
    476         Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_F','Civ1_C'};%  cell array containing the names of the fields to record
    477         Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
    478         Data.VarAttribute{1}.Role='coord_x';
    479         Data.VarAttribute{2}.Role='coord_y';
    480         Data.VarAttribute{3}.Role='vector_x';
    481         Data.VarAttribute{4}.Role='vector_y';
    482         Data.VarAttribute{5}.Role='warnflag';
    483         % case of mask
     283
     284        disp('civ1 started')
     285
     286
     287        % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
     288
     289        par_civ1.ImageA=zeros(2*SearchRange_z+1,npy,npx);
     290        par_civ1.ImageB=zeros(2*SearchRange_z+1,npy,npx);
     291
     292        for islice=ceil(par_civ1.Dz/2):par_civ1.Dz:NbSlice
     293            if par_civ1.Dz<2*SearchRange_z+1
     294            par_civ1.ImageA=circshift(par_civ1.ImageA,-par_civ1.Dz);
     295            par_civ1.ImageB=circshift(par_civ1.ImageA,-par_civ1.Dz);
     296            end
     297              for iz=1:par_civ1.Dz
     298            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield,islice));
     299            A= read_image(ImageName_A,FileType_A);
     300            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield,islice));
     301            B= read_image(ImageName_B,FileType_B);
     302            par_civ1.ImageA(2*SearchRange_z+2-iz,:,:) = A;
     303            par_civ1.ImageB(2*SearchRange_z+2-iz,:,:) = B;
     304            end
     305            % caluclate velocity data (y and v in indices, reverse to y component)
     306            [Data.Civ1_X(islice,:,:),Data.Civ1_Y(islice,:,:), utable, vtable,wtable, ctable, FF, result_conv, errormsg] = civ3D (par_civ1);
     307            if ~isempty(errormsg)
     308                disp_uvmat('ERROR',errormsg,checkrun)
     309                return
     310            end
     311
     312        end
     313        % case of mask TO ADAPT
    484314        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
    485315            if isfield(par_civ1,'NbSlice')
    … …  
    493323                par_civ1.Mask=mask; %use mask already opened
    494324            else
    495                 if exist(maskname,'file')
     325                if ~isempty(regexp(maskname,'(^http://)|(^https://)'))|| exist(maskname,'file')
    496326                    try
    497327                        par_civ1.Mask=imread(maskname);%update the mask, an store it for future use
    … …  
    510340            end
    511341        end
    512         if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume')
    513             Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
    514             Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
    515             Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];Data.Civ1_F=[];
    516             for ivol=1:NbSlice
    517                 % caluclate velocity data (y and v in indices, reverse to y component)
    518                 [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
    519                 if ~isempty(errormsg)
    520                     disp_uvmat('ERROR',errormsg,checkrun)
    521                     return
    522                 end
    523                 Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
    524                 Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
    525                 Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates
    526                 Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)];
    527                 Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)];
    528                 Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
    529                 Data.Civ1_F=[Data.Civ1_C reshape(F,[],1)];
    530             end
    531         else %usual PIV
    532             % caluclate velocity data (y and v in indices, reverse to y component)
    533             [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
    534             if ~isempty(errormsg)
    535                 disp_uvmat('ERROR',errormsg,checkrun)
    536                 return
    537             end
    538             Data.Civ1_X=reshape(xtable,[],1);
    539             Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1);
    540             Data.Civ1_U=reshape(utable,[],1);
    541             Data.Civ1_V=reshape(-vtable,[],1);
    542             Data.Civ1_C=reshape(ctable,[],1);
    543             Data.Civ1_F=reshape(F,[],1);
    544         end
    545     else% we use existing Civ1 data
    546         if exist('ncfile','var')
    547             CivFile=ncfile;
    548             [Data,tild,tild,errormsg]=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
    549             if ~isempty(errormsg)
    550                 disp_uvmat('ERROR',errormsg,checkrun)
    551                 return
    552             end
    553             [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
    554         elseif isfield(Param,'Civ1_X')
    555             Data.ListGlobalAttribute={};
    556             Data.ListVarName={};
    557             Data.VarDimName={};
    558             Data.Civ1_X=Param.Civ1_X;
    559             Data.Civ1_Y=Param.Civ1_Y;
    560             Data.Civ1_U=Param.Civ1_U;
    561             Data.Civ1_V=Param.Civ1_V;
    562             Data.Civ1_FF=Param.Civ1_FF;
    563         end
    564     end
    565    
     342        % Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
     343        % Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
     344        % Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];
     345        %
     346        %
     347        % Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
     348        % Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
     349        % Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates
     350        % Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)];
     351        % Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)];
     352        % Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
     353        % Data.Civ1_FF=[Data.Civ1_FF reshape(F,[],1)];
     354
     355    end
     356
    566357    %% Fix1
    567358    if isfield (Param.ActionInput,'Fix1')
    568         disp('fix1 started')
     359        disp('detect_false1 started')
    569360        if ~isfield (Param.ActionInput,'Civ1')% if we use existing Civ1, remove previous data beyond Civ1
    570361            Fix1_attr=find(strcmp('Fix1',Data.ListGlobalAttribute));
    … …  
    581372        end
    582373        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix1_param];
    583         Data.ListVarName=[Data.ListVarName {'Civ1_FF'}];
    584         Data.VarDimName=[Data.VarDimName {'nb_vec_1'}];
    585         nbvar=length(Data.ListVarName);
    586         Data.VarAttribute{nbvar}.Role='errorflag';
    587         Data.Civ1_FF=int8(fix(Param.ActionInput.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V));
     374        Data.Civ1_FF=uint8(detect_false(Param.ActionInput.Fix1,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V,Data.Civ1_FF));
    588375        Data.CivStage=2;
    589376    end
    … …  
    591378    if isfield (Param.ActionInput,'Patch1')
    592379        disp('patch1 started')
    593         if check_civx
    594             errormsg='Civ Matlab input needed for patch';
    595             disp_uvmat('ERROR',errormsg,checkrun)
    596             return
    597         end
    598        
     380        tstart_patch1=tic;
     381
    599382        % record the processing parameters of Patch1 as global attributes in the result nc file
    600383        list_param=fieldnames(Param.ActionInput.Patch1)';
    … …  
    606389        Data.CivStage=3;% record the new state of processing
    607390        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch1_param];
    608        
     391
    609392        % list the variables to record
    610393        nbvar=length(Data.ListVarName);
    … …  
    624407            ind_good=1:numel(Data.Civ1_X);
    625408        end
    626        
     409        if isempty(ind_good)
     410            disp_uvmat('ERROR','all vectors of civ1 are bad, check input parameters' ,checkrun)
     411            return
     412        end
     413
    627414        % perform Patch calculation using the UVMAT fct 'filter_tps'
     415
    628416        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
    629417            filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
    630418        Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
    631419        Data.Civ1_V_smooth(ind_good)=Vres;
    632         Data.Civ1_FF(ind_good)=int8(FFres);
     420        Data.Civ1_FF(ind_good)=uint8(FFres);
     421        time_patch1=toc(tstart_patch1);
    633422        disp('patch1 performed')
    634423    end
    635    
     424
    636425    %% Civ2
    637426    if isfield (Param.ActionInput,'Civ2')
    638427        disp('civ2 started')
     428        tstart_civ2=tic;
    639429        par_civ2=Param.ActionInput.Civ2;
    640         if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
    641             par_civ2.ImageA=[];
    642             par_civ2.ImageB=[];
    643             if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    644                     ImageName_A_Civ2=Param.ActionInput.RefFile;
    645                 else
     430        % read input images
     431        par_civ2.ImageA=[];
     432        par_civ2.ImageB=[];
     433        if strcmp(Param.ActionInput.ListCompareMode,'displacement')
     434            ImageName_A_Civ2=Param.ActionInput.RefFile;
     435        else
    646436            ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_civ2,[],j1_civ2);
    647             end
    648             if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield))
    649                 par_civ2.ImageA=par_civ1.ImageA;
    650             else
    651                 [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield));
    652             end
    653             ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_civ2,[],j2_civ2);
    654             if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2)
    655                 par_civ2.ImageB=par_civ1.ImageB;
    656             else
    657                 [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield));
    658             end
    659             par_civ2.ImageWidth=FileInfo_A.Width;
    660             par_civ2.ImageHeight=FileInfo_A.Height;
    661             if isfield(par_civ2,'Grid')% grid points set as input file
    662                 if ischar(par_civ2.Grid)%read the grid file if the input is a file name
    663                     par_civ2.Grid=dlmread(par_civ2.Grid);
    664                     par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
    665                 end
    666             else% automatic grid
    667                 minix=floor(par_civ2.Dx/2)-0.5;
    668                 maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
    669                 miniy=floor(par_civ2.Dy/2)-0.5;
    670                 maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
    671                 [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
    672                 par_civ2.Grid(:,1)=reshape(GridX,[],1);
    673                 par_civ2.Grid(:,2)=reshape(GridY,[],1);
    674             end
    675         end
    676        
     437        end
     438        if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield))
     439            par_civ2.ImageA=par_civ1.ImageA;
     440        else
     441            [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield));
     442        end
     443        ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_civ2,[],j2_civ2);
     444        if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2)
     445            par_civ2.ImageB=par_civ1.ImageB;
     446        else
     447            [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield));
     448        end
     449        [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A_Civ2);
     450        par_civ2.ImageWidth=FileInfo_A.Width;
     451        par_civ2.ImageHeight=FileInfo_A.Height;
     452        if isfield(par_civ2,'Grid')% grid points set as input file
     453            if ischar(par_civ2.Grid)%read the grid file if the input is a file name
     454                par_civ2.Grid=dlmread(par_civ2.Grid);
     455                par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     456            end
     457        else% automatic grid
     458            minix=floor(par_civ2.Dx/2)-0.5;
     459            maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
     460            miniy=floor(par_civ2.Dy/2)-0.5;
     461            maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
     462            [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
     463            par_civ2.Grid(:,1)=reshape(GridX,[],1);
     464            par_civ2.Grid(:,2)=reshape(GridY,[],1);
     465        end
     466        % end
     467
    677468        % get the guess from patch1 or patch2 (case 'CheckCiv3')
    678         if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
    679             if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
    680                 SubRange= Data.Civ2_SubRange;
    681                 NbCentres=Data.Civ2_NbCentres;
    682                 Coord_tps=Data.Civ2_Coord_tps;
    683                 U_tps=Data.Civ2_U_tps;
    684                 V_tps=Data.Civ2_V_tps;
    685                 CivStage=Data.CivStage;%store the current CivStage
    686                 Civ1_Dt=Data.Civ2_Dt;
    687                 Data=[];%reinitialise the result structure Data
    688                 Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
    689                 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
    690                 Data.Program='civ_series';
    691                 Data.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data
    692                 Data.ListVarName={};
    693                 Data.VarDimName={};
    694             else % get the guess from patch1
    695                 SubRange= Data.Civ1_SubRange;
    696                 NbCentres=Data.Civ1_NbCentres;
    697                 Coord_tps=Data.Civ1_Coord_tps;
    698                 U_tps=Data.Civ1_U_tps;
    699                 V_tps=Data.Civ1_V_tps;
    700                 Civ1_Dt=Data.Civ1_Dt;
    701                 Data.CivStage=4;
    702             end
    703         else
    704             SubRange= par_civ2.Civ1_SubRange;
    705             NbCentres=par_civ2.Civ1_NbCentres;
    706             Coord_tps=par_civ2.Civ1_Coord_tps;
    707             U_tps=par_civ2.Civ1_U_tps;
    708             V_tps=par_civ2.Civ1_V_tps;
    709             Civ1_Dt=par_civ2.Civ1_Dt;
    710             Civ2_Dt=par_civ2.Civ1_Dt;
     469
     470        if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
     471            SubRange= Data.Civ2_SubRange;
     472            NbCentres=Data.Civ2_NbCentres;
     473            Coord_tps=Data.Civ2_Coord_tps;
     474            U_tps=Data.Civ2_U_tps;
     475            V_tps=Data.Civ2_V_tps;
     476            CivStage=Data.CivStage;%store the current CivStage
     477            Civ1_Dt=Data.Civ2_Dt;
     478            Data=[];%reinitialise the result structure Data
     479            Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
     480            Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
     481            Data.Program='civ_series';
     482            Data.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data
    711483            Data.ListVarName={};
    712484            Data.VarDimName={};
    713         end
     485        else % get the guess from patch1
     486            SubRange= Data.Civ1_SubRange;
     487            NbCentres=Data.Civ1_NbCentres;
     488            Coord_tps=Data.Civ1_Coord_tps;
     489            U_tps=Data.Civ1_U_tps;
     490            V_tps=Data.Civ1_V_tps;
     491            Civ1_Dt=Data.Civ1_Dt;
     492            Data.CivStage=4;
     493        end
     494        % else
     495        %     SubRange= par_civ2.Civ1_SubRange;
     496        %     NbCentres=par_civ2.Civ1_NbCentres;
     497        %     Coord_tps=par_civ2.Civ1_Coord_tps;
     498        %     U_tps=par_civ2.Civ1_U_tps;
     499        %     V_tps=par_civ2.Civ1_V_tps;
     500        %     Civ1_Dt=par_civ2.Civ1_Dt;
     501        %     Civ2_Dt=par_civ2.Civ1_Dt;
     502        %     Data.ListVarName={};
     503        %     Data.VarDimName={};
     504        % end
    714505        Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
    715506        Shifty=zeros(size(par_civ2.Grid,1),1);
    … …  
    742533            end
    743534        end
    744         if par_civ2.CheckMask&&~isempty(par_civ2.Mask)       
     535        if par_civ2.CheckMask&&~isempty(par_civ2.Mask)
    745536            if isfield(par_civ2,'NbSlice')
    746537                [RootPath_mask,SubDir_mask,RootFile_mask,i1_mask,i2_mask,j1_mask,j2_mask,Ext_mask]=fileparts_uvmat(Param.ActionInput.Civ2.Mask);
    … …  
    754545            else
    755546                if exist(maskname,'file')
    756                 try
    757                     par_civ2.Mask=imread(maskname);%update the mask, an store it for future use
    758                 catch ME
    759                     if ~isempty(ME.message)
    760                         errormsg=['error reading input image: ' ME.message];
    761                         disp_uvmat('ERROR',errormsg,checkrun)
    762                         return
     547                    try
     548                        par_civ2.Mask=imread(maskname);%update the mask, an store it for future use
     549                    catch ME
     550                        if ~isempty(ME.message)
     551                            errormsg=['error reading input image: ' ME.message];
     552                            disp_uvmat('ERROR',errormsg,checkrun)
     553                            return
     554                        end
    763555                    end
    764                 end
    765556                else
    766557                    par_civ2.Mask=[];
    … …  
    771562        end
    772563
    773         if CheckInputFile % else Dt given by par_civ2
    774             if strcmp(Param.ActionInput.ListCompareMode,'displacement')
    775                 Civ1_Dt=1;
    776                 Civ2_Dt=1;
    777             else
    778                 Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1);
    779             end
    780         end
     564
     565        if strcmp(Param.ActionInput.ListCompareMode,'displacement')
     566            Civ1_Dt=1;
     567            Civ2_Dt=1;
     568        else
     569            Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1);
     570        end
     571
    781572        par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
    782573        % shift the grid points by half the expected shift to provide the correlation box position in image A
    … …  
    788579            par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
    789580        end
    790        
     581
    791582        % calculate velocity data (y and v in image indices, reverse to y component)
     583
    792584        [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2);
    793        
     585
    794586        list_param=(fieldnames(Param.ActionInput.Civ2))';
    795587        list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list
    … …  
    812604        end
    813605        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param];
    814        
     606
    815607        nbvar=numel(Data.ListVarName);
    816608        % define the Civ2 variable (if Civ2 data are not replaced from previous calculation)
    817609        if isempty(find(strcmp('Civ2_X',Data.ListVarName),1))
    818             Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];%  cell array containing the names of the fields to record
     610            Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_C','Civ2_FF'}];%  cell array containing the names of the fields to record
    819611            Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
    820612            Data.VarAttribute{nbvar+1}.Role='coord_x';
    … …  
    822614            Data.VarAttribute{nbvar+3}.Role='vector_x';
    823615            Data.VarAttribute{nbvar+4}.Role='vector_y';
    824             Data.VarAttribute{nbvar+5}.Role='warnflag';
     616            Data.VarAttribute{nbvar+5}.Role='ancillary';
     617            Data.VarAttribute{nbvar+6}.Role='errorflag';
    825618        end
    826619        Data.Civ2_X=reshape(xtable,[],1);
    … …  
    829622        Data.Civ2_V=reshape(-vtable,[],1);
    830623        Data.Civ2_C=reshape(ctable,[],1);
    831         Data.Civ2_F=reshape(F,[],1);
     624        Data.Civ2_FF=reshape(F,[],1);
    832625        disp('civ2 performed')
     626        time_civ2=toc(tstart_civ2);
    833627    elseif ~isfield(Data,'ListVarName') % we start there, using existing Civ2 data
    834628        if exist('ncfile','var')
    835629            CivFile=ncfile;
    836             [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
     630            [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and detect_false1 data in the existing netcdf file
    837631            if ~isempty(errormsg)
    838632                disp_uvmat('ERROR',errormsg,checkrun)
    839633                return
    840634            end
    841         elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)
    842             Data.ListGlobalAttribute={};
    843             Data.ListVarName={};
    844             Data.VarDimName={};
    845             Data.Civ2_X=Param.Civ2_X;
    846             Data.Civ2_Y=Param.Civ2_Y;
    847             Data.Civ2_U=Param.Civ2_U;
    848             Data.Civ2_V=Param.Civ2_V;
    849             Data.Civ2_FF=Param.Civ2_FF;
    850         end
    851     end
    852    
     635            %         elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)
     636            %             Data.ListGlobalAttribute={};
     637            %             Data.ListVarName={};
     638            %             Data.VarDimName={};
     639            %             Data.Civ2_X=Param.Civ2_X;
     640            %             Data.Civ2_Y=Param.Civ2_Y;
     641            %             Data.Civ2_U=Param.Civ2_U;
     642            %             Data.Civ2_V=Param.Civ2_V;
     643            %             Data.Civ2_FF=Param.Civ2_FF;
     644        end
     645    end
     646
    853647    %% Fix2
    854648    if isfield (Param.ActionInput,'Fix2')
    855         disp('fix2 started')
     649        disp('detect_false2 started')
    856650        list_param=fieldnames(Param.ActionInput.Fix2)';
    857651        Fix2_param=regexprep(list_param,'^.+','Fix2_$0');% insert 'Fix1_' before  each string in ListFixParam
    … …  
    861655        end
    862656        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];
    863         %
    864         %         ListFixParam=fieldnames(Param.ActionInput.Fix2);
    865         %         for ilist=1:length(ListFixParam)
    866         %             ParamName=ListFixParam{ilist};
    867         %             ListName=['Fix2_' ParamName];
    868         %             eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
    869         %             eval(['Data.' ListName '=Param.ActionInput.Fix2.' ParamName ';'])
    870         %         end
    871         if check_civx
    872             if ~isfield(Data,'fix2')
    873                 Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2'];
    874                 Data.fix2=1;
    875                 Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}];
    876                 Data.VarDimName=[Data.VarDimName {'nb_vectors2'}];
    877             end
    878             Data.vec_FixFlag=fix(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y);
    879         else
    880             Data.ListVarName=[Data.ListVarName {'Civ2_FF'}];
    881             Data.VarDimName=[Data.VarDimName {'nb_vec_2'}];
    882             nbvar=length(Data.ListVarName);
    883             Data.VarAttribute{nbvar}.Role='errorflag';
    884             Data.Civ2_FF=double(fix(Param.ActionInput.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V));
    885             Data.CivStage=Data.CivStage+1;
    886         end
    887     end
    888    
     657        Data.Civ2_FF=double(detect_false(Param.ActionInput.Fix2,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V,Data.Civ2_FF));
     658        Data.CivStage=Data.CivStage+1;
     659    end
     660
    889661    %% Patch2
    890662    if isfield (Param.ActionInput,'Patch2')
    891663        disp('patch2 started')
     664        tstart_patch2=tic;
    892665        list_param=fieldnames(Param.ActionInput.Patch2)';
    893666        list_param(strcmp('TestPatch2',list_param))=[];% remove the parameter TestCiv1 from the list
    … …  
    898671        end
    899672        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch2_param];
    900        
     673
    901674        nbvar=length(Data.ListVarName);
    902675        Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
    903676        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
    904677            {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
    905        
     678
    906679        Data.VarAttribute{nbvar+1}.Role='vector_x';
    907680        Data.VarAttribute{nbvar+2}.Role='vector_y';
    … …  
    916689            ind_good=1:numel(Data.Civ2_X);
    917690        end
     691        if isempty(ind_good)
     692            disp_uvmat('ERROR','all vectors of civ2 are bad, check input parameters' ,checkrun)
     693            return
     694        end
     695
    918696        [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=...
    919697            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);
    … …  
    922700        Data.Civ2_FF(ind_good)=FFres;
    923701        Data.CivStage=Data.CivStage+1;
     702        time_patch2=toc(tstart_patch2);
    924703        disp('patch2 performed')
    925704    end
    926    
     705
    927706    %% write result in a netcdf file if requested
    928     if CheckOutputFile
    929         errormsg=struct2nc(ncfile_out,Data);
    930         if isempty(errormsg)
    931             disp([ncfile_out ' written'])
    932             %[success,msg] = fileattrib(ncfile_out ,'+w','g');% done in struct2nc
    933         else
    934             disp(errormsg)
    935         end
    936         disp(['ellapsed time ' num2str(toc/60,2) ' minutes'])
    937     end
     707    % if CheckOutputFile
     708    errormsg=struct2nc(ncfile_out,Data);
     709    if isempty(errormsg)
     710        disp([ncfile_out ' written'])
     711        %[success,msg] = fileattrib(ncfile_out ,'+w','g');% done in struct2nc
     712    else
     713        disp(errormsg)
     714    end
     715    time_total=toc(tstart);
     716    disp(['ellapsed time ' num2str(time_total/60,2) ' minutes'])
     717    disp(['time image reading ' num2str(time_input,2) ' s'])
     718    disp(['time civ1 ' num2str(time_civ1,2) ' s'])
     719    disp(['time patch1 ' num2str(time_patch1,2) ' s'])
     720    disp(['time civ2 ' num2str(time_civ2,2) ' s'])
     721    disp(['time patch2 ' num2str(time_patch2,2) ' s'])
     722    disp(['time other ' num2str((time_total-time_input-time_civ1-time_patch1-time_civ2-time_patch2),2) ' s'])
     723    % end
    938724end
    939725
    … …  
    971757%  .DVDY
    972758
    973 function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ)
     759function [xtable,ytable,utable,vtable,wtable,ctable,FF,result_conv,errormsg] = civ3D (par_civ)
    974760
    975761%% prepare measurement grid
    976 if isfield(par_civ,'Grid')% grid points set as input, central positions of the sub-images in image A
    977     if ischar(par_civ.Grid)%read the grid file if the input is a file name (grid in x, y image coordinates)
    978         par_civ.Grid=dlmread(par_civ.Grid);
    979         par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
    980     end
    981     % else par_civ.Grid is already an array, no action here
    982 else% automatic grid in x, y image coordinates
    983     minix=floor(par_civ.Dx/2)-0.5;
    984     maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
    985     miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
    986     maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
    987     [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
    988     par_civ.Grid(:,1)=reshape(GridX,[],1);
    989     par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
    990 end
    991 nbvec=size(par_civ.Grid,1);
     762
     763minix=floor(par_civ.Dx/2)-0.5;
     764maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
     765miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
     766maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
     767[GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     768par_civ.Grid(:,:,1)=GridX;
     769par_civ.Grid(:,:,2)=GridY;% increases with array index,
     770[nbvec_y,nbvec_x,~]=size(par_civ.Grid);
     771%
     772%
     773% minix=floor(par_civ.Dx/2)-0.5;
     774%     maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
     775%     miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
     776%     maxiy=miniy+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
     777%     [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     778%     par_civ.Grid(:,1)=reshape(GridX,[],1);
     779%     par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
    992780
    993781%% prepare correlation and search boxes
    … …  
    996784isx2=floor(par_civ.SearchBoxSize(1)/2);
    997785isy2=floor(par_civ.SearchBoxSize(2)/2);
     786isz2=floor(par_civ.SearchBoxSize(3)/2);
     787kref=isz2+1;%middle index of the z slice
    998788shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value
    999789shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
    1000790if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
    1001     shiftx=shiftx*ones(nbvec,1);
    1002     shifty=shifty*ones(nbvec,1);
     791    shiftx=shiftx*ones(nbvec_y,nbvec_x,1);
     792    shifty=shifty*ones(nbvec_y,nbvec_x,1);
    1003793end
    1004794
    1005795%% Array initialisation and default output  if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation)
    1006 xtable=round(par_civ.Grid(:,1)+0.5)-0.5;
    1007 ytable=round(par_civ.ImageHeight-par_civ.Grid(:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
     796xtable=round(par_civ.Grid(:,:,1)+0.5)-0.5;
     797ytable=round(par_civ.ImageHeight-par_civ.Grid(:,:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
    1008798utable=shiftx;%zeros(nbvec,1);
    1009799vtable=shifty;%zeros(nbvec,1);
    1010 ctable=zeros(nbvec,1);
    1011 F=zeros(nbvec,1);
     800wtable=zeros(size(utable));
     801ctable=zeros(nbvec_y,nbvec_x,1);
     802FF=zeros(nbvec_y,nbvec_x,1);
    1012803result_conv=[];
    1013804errormsg='';
    1014805
    1015806%% prepare mask
    1016 if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
    1017     if strcmp(par_civ.Mask,'all')
    1018         return    % get the grid only, no civ calculation
    1019     elseif ischar(par_civ.Mask)
    1020         par_civ.Mask=imread(par_civ.Mask);% read the mask if not allready done
    1021     end
    1022 end
    1023807check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
    1024808check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
    1025809
    1026 par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
    1027 par_civ.ImageB=sum(double(par_civ.ImageB),3);
    1028 [npy_ima npx_ima]=size(par_civ.ImageA);
    1029 if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima])
     810[npz,npy_ima,npx_ima]=size(par_civ.ImageA);
     811if ~isequal(size(par_civ.ImageB),[npz npy_ima npx_ima])
    1030812    errormsg='image pair with unequal size';
    1031813    return
    … …  
    1033815
    1034816%% Apply mask
    1035 % Convention for mask IDEAS TO IMPLEMENT ?
     817% Convention for mask, IDEAS NOT IMPLEMENTED
    1036818% mask >200 : velocity calculated
    1037819%  200 >=mask>150;velocity not calculated, interpolation allowed (bad spots)
    … …  
    1040822%  20>=mask: velocity=0
    1041823checkmask=0;
    1042 MinA=min(min(par_civ.ImageA));
     824MinA=min(min(min(par_civ.ImageA)));
    1043825%MinB=min(min(par_civ.ImageB));
    1044826%check_undefined=false(size(par_civ.ImageA));
    … …  
    1050832    end
    1051833    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
    1052     %     par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined  area
    1053     %     par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined  area
    1054834end
    1055835
    … …  
    1065845
    1066846if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given)
    1067     for ivec=1:nbvec
    1068         iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
    1069         jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);%  j index  for the middle of the correlation box in the image A
    1070         F(ivec)=0;
    1071         subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
    1072         subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
    1073         subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
    1074         subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
    1075         image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
    1076         image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
    1077         check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
    1078         check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
    1079         check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
    1080         check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
    1081         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
    1082         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
    1083         if checkmask
    1084             mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
    1085             mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
    1086             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
    1087             mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B
    1088             sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
    1089             if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
    1090                 F(ivec)=3; %
    1091                 utable(ivec)=0;
    1092                 vtable(ivec)=0;
     847    for ivec_x=1:nbvec_x
     848        for ivec_y=1:nbvec_y
     849            ivec_y
     850            iref=round(par_civ.Grid(ivec_y,ivec_x,1)+0.5)% xindex on the image A for the middle of the correlation box
     851            jref=round(par_civ.ImageHeight-par_civ.Grid(ivec_y,ivec_x,2)+0.5)%  j index  for the middle of the correlation box in the image A
     852            subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
     853            subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
     854            subrange2_x=iref+shiftx(ivec_y,ivec_x)-isx2:iref+shiftx(ivec_y,ivec_x)+isx2;%x indices defining the second subimage
     855            subrange2_y=jref+shifty(ivec_y,ivec_x)-isy2:jref+shifty(ivec_y,ivec_x)+isy2;%y indices defining the second subimage
     856            image1_crop=MinA*ones(npz,numel(subrange1_y),numel(subrange1_x));% default value=min of image A
     857            image2_crop=MinA*ones(npz,numel(subrange2_y),numel(subrange2_x));% default value=min of image A
     858            check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
     859            check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     860            check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
     861            check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
     862            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
     863            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
     864            if checkmask
     865                mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
     866                mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
     867                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
     868                mask2_crop(check2_y,check2_x)=check_undefined(subrange2_y(check2_y),subrange2_x(check2_x));%extract a mask subimage (search box) from image B
     869                sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
     870                if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
     871                    FF(ivec_y,ivec_x)=1; %
     872                    utable(ivec_y,ivec_x)=NaN;
     873                    vtable(ivec_y,ivec_x)=NaN;
     874                else
     875                    FF(ivec_y,ivec_x)=0;
     876                    image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
     877                    image2_crop=image2_crop.*~mask2_crop;
     878                    image1_mean=mean(mean(image1_crop))/(1-sizemask);
     879                    image2_mean=mean(mean(image2_crop))/(1-sizemask);
     880                end
    1093881            else
    1094                 image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
    1095                 image2_crop=image2_crop.*~mask2_crop;
    1096                 image1_mean=mean(mean(image1_crop))/(1-sizemask);
    1097                 image2_mean=mean(mean(image2_crop))/(1-sizemask);
    1098             end
    1099         else
    1100             image1_mean=mean(mean(image1_crop));
    1101             image2_mean=mean(mean(image2_crop));
    1102         end
    1103         %threshold on image minimum
    1104         if F(ivec)~=3
    1105             if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
    1106                 F(ivec)=3;
    1107                 %threshold on image maximum
    1108             elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
    1109                 F(ivec)=3;
    1110             end
    1111             if F(ivec)==3
    1112                 utable(ivec)=0;
    1113                 vtable(ivec)=0;
    1114             else
    1115                 %mask
    1116                 if checkmask
    1117                     image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
    1118                     image2_crop=(image2_crop-image2_mean).*~mask2_crop;
     882                image1_mean=mean(mean(image1_crop));
     883                image2_mean=mean(mean(image2_crop));
     884            end
     885            %threshold on image minimum
     886            if FF(ivec_y,ivec_x)==0
     887                if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
     888                    FF(ivec_y,ivec_x)=1;
     889                    %threshold on image maximum
     890                elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
     891                    FF(ivec_y,ivec_x)=1;
     892                end
     893                if FF(ivec_y,ivec_x)==1
     894                    utable(ivec_y,ivec_x)=NaN;
     895                    vtable(ivec_y,ivec_x)=NaN;
    1119896                else
    1120                     image1_crop=(image1_crop-image1_mean);
    1121                     image2_crop=(image2_crop-image2_mean);
    1122                 end
    1123                 %deformation
    1124                 if CheckDeformation
    1125                     xi=(1:mesh:size(image1_crop,2));
    1126                     yi=(1:mesh:size(image1_crop,1))';
    1127                     [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
    1128                     XIant=XI-par_civ.DUDX(ivec)*XI+par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
    1129                     YIant=YI+par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
    1130                     image1_crop=interp2(image1_crop,XIant,YIant);
    1131                     image1_crop(isnan(image1_crop))=0;
    1132                     xi=(1:mesh:size(image2_crop,2));
    1133                     yi=(1:mesh:size(image2_crop,1))';
    1134                     image2_crop=interp2(image2_crop,xi,yi,'*spline');
    1135                     image2_crop(isnan(image2_crop))=0;
    1136                 end
    1137                 sum_square=sum(sum(image1_crop.*image1_crop));
    1138                 %reference: Oliver Pust, PIV: Direct Cross-Correlation
    1139                 result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
    1140                 corrmax= max(max(result_conv));
    1141                 result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
    1142                 %Find the correlation max, at 255
    1143                 [y,x] = find(result_conv==255,1);
    1144                 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
    1145                 sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
    1146                 sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
    1147                 if ~isempty(y) && ~isempty(x)
    1148                     try
    1149                         if par_civ.CorrSmooth==1
    1150                             [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
    1151                         elseif par_civ.CorrSmooth==2
    1152                             [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
    1153                         else
    1154                             [vector,F(ivec)] = quadr_fit(result_conv,x,y);
    1155                         end
    1156                         utable(ivec)=vector(1)*mesh+shiftx(ivec);
    1157                         vtable(ivec)=vector(2)*mesh+shifty(ivec);
    1158                         xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
    1159                         ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
    1160                         iref=round(xtable(ivec)+0.5);% nearest image index for the middle of the vector
    1161                         jref=round(ytable(ivec)+0.5);
    1162                         % eliminate vectors located in the mask
    1163                         if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
    1164                             utable(ivec)=0;
    1165                             vtable(ivec)=0;
    1166                             F(ivec)=3;
    1167                         end
    1168                         ctable(ivec)=corrmax/sum_square;% correlation value
    1169                     catch ME
    1170                         F(ivec)=3;
    1171                         disp(ME.message)
     897                    %mask
     898                    if checkmask
     899                        image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
     900                        image2_crop=(image2_crop-image2_mean).*~mask2_crop;
     901                    else
     902                        image1_crop=(image1_crop-image1_mean);
     903                        image2_crop=(image2_crop-image2_mean);
    1172904                    end
    1173                 else
    1174                     F(ivec)=3;
     905
     906               
     907                    %reference: Oliver Pust, PIV: Direct Cross-Correlation
     908                    for kz=1:par_civ.SearchBoxSize(3)
     909                        subima2=squeeze(image2_crop(kz,:,:));
     910                        subima1=squeeze(image1_crop(kref,:,:));
     911                        correl_xy=conv2(subima2,flip(flip(subima1,2),1),'valid');
     912                          result_conv(kz,:,:)= correl_xy;
     913                        max_xy(kz)=max(max(correl_xy));
     914                    [xk(kz),yk(kz)]=find(correl_xy==max_xy(kz),1);
     915               
     916                    end
     917                    [corrmax,z]=max(max_xy);
     918               
     919                    x=xk(z);
     920                    y=yk(z);
     921                    result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
     922                    subimage2_crop=squeeze(image2_crop(z,y:y+2*iby2/mesh,x:x+2*ibx2/mesh));%subimage of image 2 corresponding to the optimum displacement of first image
     923                    sum_square=sum(sum(squeeze(image1_crop(z,:,:).*image1_crop(z,:,:))));
     924                    sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
     925                    sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
     926                    if ~isempty(y) && ~isempty(x)
     927       
     928                            if par_civ.CorrSmooth==1
     929                                [vector,FF(ivec_y,ivec_x)] = SUBPIXGAUSS (result_conv(z,:,:),x,y);%TODO: improve by max optimisation along z
     930                            elseif par_civ.CorrSmooth==2
     931                                [vector,FF(ivec_y,ivec_x)] = SUBPIX2DGAUSS (result_conv(z,:,:),x,y);
     932                            else
     933                                [vector,FF(ivec_y,ivec_x)] = quadr_fit(result_conv(z,:,:),x,y);
     934                            end
     935                            utable(ivec_y,ivec_x)=vector(1)*mesh+shiftx(ivec_y,ivec_x);
     936                            vtable(ivec_y,ivec_x)=vector(2)*mesh+shifty(ivec_y,ivec_x);
     937                            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)
     938                            ytable(ivec_y,ivec_x)=jref+vtable(ivec_y,ivec_x)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     939                            iref=round(xtable(ivec_y,ivec_x)+0.5);% nearest image index for the middle of the vector
     940                            jref=round(ytable(ivec_y,ivec_x)+0.5);
     941                            wtable(ivec_y,ivec_x)=z-kref;
     942                            % eliminate vectors located in the mask
     943                            if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
     944                                utable(ivec_y,ivec_x)=0;
     945                                vtable(ivec_y,ivec_x)=0;
     946                                FF(ivec_y,ivec_x)=1;
     947                            end
     948                            ctable(ivec_y,ivec_x)=corrmax/sum_square;% correlation value
     949
     950                    else
     951                        FF(ivec_y,ivec_x)=1;
     952                    end
    1175953                end
    1176954            end
    … …  
    1204982    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    1205983else
    1206     F=-2; % warning flag for vector truncated by the limited search box
     984    F=1; % warning flag for vector truncated by the limited search box
    1207985end
    1208986peakx=x;
    … …  
    1213991    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    1214992else
    1215     F=-2; % warning flag for vector truncated by the limited search box
     993    F=1; % warning flag for vector truncated by the limited search box
    1216994end
    1217995vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
    … …  
    12221000%------------------------------------------------------------------------
    12231001% vector=[0 0]; %default
    1224 F=-2;
     1002F=1;
    12251003peaky=y;
    12261004peakx=x;
    … …  
    12651043[npy,npx]=size(result_conv);
    12661044if x<4 || y<4 || npx-x<4 ||npy-y <4
    1267     F=-2;
     1045    F=1;
    12681046    vector=[x y];
    12691047else
    … …  
    12981076end
    12991077
    1300 %'RUN_FIX': function for fixing velocity fields:
    1301 %-----------------------------------------------
    1302 % RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)
    1303 %
    1304 %filename: name of the netcdf file (used as input and output)
    1305 %field: structure specifying the names of the fields to fix (depending on civ1 or civ2)
    1306 %.vel_type='civ1' or 'civ2';
    1307 %.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1);
    1308 %.fixflag=name of fix flag variable ('vec_FixFlag' for civ1)
    1309 %flagindex: flag specifying which values of vec_f are removed:
    1310 % if flagindex(1)=1: vec_f=-2 vectors are removed
    1311 % if flagindex(2)=1: vec_f=3 vectors are removed
    1312 % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
    1313 %iter=1 for civ1 fields and iter=2 for civ2 fields
    1314 %thresh_vecC: threshold in the image correlation vec_C
    1315 %flag_mask: =1 mask used to remove vectors (0 else)
    1316 %maskname: name of the mask image file for fix
    1317 %thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
    1318 %inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
    1319 %fileref: .nc file name for a reference velocity (='': refrence 0 used)
    1320 %fieldref: 'civ1','filter1'...feld used in fileref
    1321 
    1322 function FF=fix(Param,F,C,U,V,X,Y)
    1323 FF=zeros(size(F));%default
    1324 
    1325 %criterium on warn flags
    1326 FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'};
    1327 FlagVal=[-2 2 3 4];
    1328 for iflag=1:numel(FlagName)
    1329     if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag})
    1330         FF=(FF==1| F==FlagVal(iflag));
    1331     end
    1332 end
    1333 %criterium on correlation values
     1078
     1079function FF=detect_false(Param,C,U,V,FFIn)
     1080FF=FFIn;%default, good vectors
     1081% FF=1, for correlation max at edge, not set in this function
     1082% FF=2, for too small correlation
     1083% FF=3, for velocity outside bounds
     1084% FF=4 for exclusion by difference with the smoothed field, not set in this function
     1085
    13341086if isfield (Param,'MinCorr')
    1335     FF=FF==1 | C<Param.MinCorr;
     1087     FF(C<Param.MinCorr & FFIn==0)=2;
    13361088end
    13371089if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
    13381090    Umod= U.*U+V.*V;
    13391091    if isfield (Param,'MinVel')&&~isempty(Param.MinVel)
    1340         FF=FF==1 | Umod<(Param.MinVel*Param.MinVel);
     1092        U2Min=Param.MinVel*Param.MinVel;
     1093        FF(Umod<U2Min & FFIn==0)=3;
    13411094    end
    13421095    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
    1343         FF=FF==1 | Umod>(Param.MaxVel*Param.MaxVel);
    1344     end
    1345 end
    1346 
     1096         U2Max=Param.MaxVel*Param.MaxVel;
     1097        FF(Umod>U2Max & FFIn==0)=3;
     1098    end
     1099end
    13471100
    13481101%------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.