source: trunk/src/series/civ_series.m @ 995

Last change on this file since 995 was 995, checked in by sommeria, 8 years ago

civ_series

File size: 67.1 KB
Line 
1%'civ_series': 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
6%------------------------------------------------------------------------
7% function [Data,errormsg,result_conv]= civ_series(Param)
8%
9%OUTPUT
10% Data=structure containing the PIV results and information on the processing parameters
11% errormsg=error message char string, decd ..fault=''
12% resul_conv: image inter-correlation function for the last grid point (used for tests)
13%
14%INPUT:
15% Param: Matlab structure of input  parameters
16%     Param contains info of the GUI series using the fct read_GUI.
17%     Param.Action.RUN = 0 (to set the status of the GUI series) or =1 to RUN the computation
18%     Param.InputTable: sets the input file(s)
19%           if absent, the fct looks for input data in Param.ActionInput     (test mode)
20%     Param.OutputSubDir: sets the folder name of output file(s,
21%           if absent no file is produced, result in the output structure Data (test mode)
22%     Param.ActionInput: substructure with the parameters provided by the GUI civ_input
23%                      .Civ1: parameters for civ1
24%                      .Fix1: parameters for fix1
25%                      .Patch1:
26%                      .Civ2: for civ2
27%                      .Fix2:
28%                      .Patch2:
29
30%=======================================================================
31% Copyright 2008-2017, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
32%   http://www.legi.grenoble-inp.fr
33%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
34%
35%     This file is part of the toolbox UVMAT.
36%
37%     UVMAT is free software; you can redistribute it and/or modify
38%     it under the terms of the GNU General Public License as published
39%     by the Free Software Foundation; either version 2 of the license,
40%     or (at your option) any later version.
41%
42%     UVMAT is distributed in the hope that it will be useful,
43%     but WITHOUT ANY WARRANTY; without even the implied warranty of
44%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
45%     GNU General Public License (see LICENSE.txt) for more details.
46%=======================================================================
47
48function [Data,errormsg,result_conv]= civ_series(Param)
49errormsg='';
50
51%% set the input elements needed on the GUI series when the action is selected in the menu ActionName or InputTable refreshed
52if isstruct(Param) && isequal(Param.Action.RUN,0)% function activated from the GUI series but not RUN
53    path_series=fileparts(which('series'));
54    addpath(fullfile(path_series,'series'))
55    Data=civ_input(Param);% introduce the civ parameters using the GUI civ_input
56    if isempty(Data)
57        Data=Param;% if  civ_input has been cancelled, keep previous parameters
58    end
59    Data.Program=mfilename;%gives the name of the current function
60    Data.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
61    Data.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
62    %     if isfield(Data,'ActionInput') && isfield(Data.ActionInput,'PairIndices')&& strcmp(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
63    %         Data.Desable_j_index='on';% hide the j index in series (set by the pair choice in civ_input)
64    %     end
65    Data.NbSlice='off'; %nbre of slices ('off' by default)
66    Data.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
67    Data.FieldName='on';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
68    Data.FieldTransform = 'off';%can use a transform function
69    Data.ProjObject='off';%can use projection object(option 'off'/'on',
70    Data.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
71    Data.OutputDirExt='.civ';%set the output dir extension
72    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)
73    Data.OutputFileMode='NbInput_i';% one output file expected per value of i index (used for waitbar)
74    Data.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1)
75    if isequal(Data.ActionInput.PairIndices.ListPairMode,'pair j1-j2')
76        if isfield(Data.ActionInput.PairIndices,'ListPairCiv2')
77            str_civ=Data.ActionInput.PairIndices.ListPairCiv2;
78        else
79            str_civ=Data.ActionInput.PairIndices.ListPairCiv1;
80        end
81        r=regexp(str_civ,'^j= (?<num1>[a-z])-(?<num2>[a-z])','names');
82        if isempty(r)
83            r=regexp(str_civ,'^j= (?<num1>[A-Z])-(?<num2>[A-Z])','names');
84            if isempty(r)
85                r=regexp(str_civ,'^j= (?<num1>\d+)-(?<num2>\d+)','names');
86            end
87        end
88        if ~isempty(r)
89            Data.j_index_1=stra2num(r.num1);
90            Data.j_index_2=stra2num(r.num2);
91        end
92    end
93    return
94end
95
96%% read input parameters from an xml file if input is a file name (batch mode)
97checkrun=1;
98if ischar(Param)
99    Param=xml2struct(Param);% read Param as input file (batch case)
100    checkrun=0;
101end
102
103%% test input
104if ~isfield(Param,'ActionInput')
105    disp_uvmat('ERROR','no parameter set for PIV',checkrun)
106    return
107end
108iview_A=0;%default values
109NbField=1;
110RUNHandle=[];
111CheckInputFile=isfield(Param,'InputTable');%= 1 in test use for TestCiv (no nc file involved)
112CheckOutputFile=isfield(Param,'OutputSubDir');%= 1 in test use for TestPatch (no nc file produced)
113
114%% input files and indexing (skipped in Test mode)
115if CheckInputFile
116    hseries=findobj(allchild(0),'Tag','series');
117    RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
118    WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
119    MaxIndex_i=Param.IndexRange.MaxIndex_i;
120    MinIndex_i=Param.IndexRange.MinIndex_i;
121    MaxIndex_j=ones(size(MaxIndex_i));MinIndex_j=ones(size(MinIndex_i));
122    if isfield(Param.IndexRange,'MaxIndex_j')&& isfield(Param.IndexRange,'MinIndex_j')
123        MaxIndex_j=Param.IndexRange.MaxIndex_j;
124        MinIndex_j=Param.IndexRange.MinIndex_j;
125    end
126    if isfield(Param,'InputTable')
127        [tild,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
128        iview_A=0;% series index (iview) for the first image series
129        iview_B=0;% series index (iview) for the second image series (only non zero for option 'shift' comparing two image series )
130        if Param.ActionInput.CheckCiv1
131            iview_A=1;% usual PIV, the image series is on the first line of the table
132        elseif Param.ActionInput.CheckCiv2 % civ2 is performed without Civ1, a netcdf file series is needed in the first table line
133            iview_A=2;% the second line is used for the input images of Civ2
134        end
135        if strcmp(Param.ActionInput.ListCompareMode,'shift')
136            iview_B=iview_A+1; % the second image series is on the next line of the input table
137        end
138        if iview_A~=0
139            RootPath_A=Param.InputTable{iview_A,1};
140            RootFile_A=Param.InputTable{iview_A,3};
141            SubDir_A=Param.InputTable{iview_A,2};
142            NomType_A=Param.InputTable{iview_A,4};
143            FileExt_A=Param.InputTable{iview_A,5};
144            if iview_B==0
145                iview_B=iview_A;% the second image series is the same as the first
146            end
147            RootPath_B=Param.InputTable{iview_B,1};
148            RootFile_B=Param.InputTable{iview_B,3};
149            SubDir_B=Param.InputTable{iview_B,2};
150            NomType_B=Param.InputTable{iview_B,4};
151            FileExt_B=Param.InputTable{iview_B,5};
152        end
153       
154        PairCiv2='';
155        switch Param.ActionInput.ListCompareMode
156            case 'PIV'
157                PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
158                if isfield(Param.ActionInput.PairIndices,'ListPairCiv2')
159                    PairCiv2=Param.ActionInput.PairIndices.ListPairCiv2;%string which determines the civ2 pair
160                end
161                if iview_A==1% if Civ1 is performed
162                    [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
163                        find_pair_indices(PairCiv1,i1_series{1},j1_series{1},MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j);
164                    if ~isempty(PairCiv2)
165                        [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds_Civ2]=...
166                            find_pair_indices(PairCiv2,i1_series{1},j1_series{1},MinIndex_i(1),MaxIndex_i(1),MinIndex_j(1),MaxIndex_j(1));
167                        check_bounds=check_bounds | check_bounds_Civ2;
168                    end
169                else% we start from an existing Civ1 file
170                    i1_series_Civ1=i1_series{1};
171                    i2_series_Civ1=i2_series{1};
172                    j1_series_Civ1=j1_series{1};
173                    j2_series_Civ1=j2_series{1};
174                    NomTypeNc=Param.InputTable{1,4};
175                    if ~isempty(PairCiv2)
176                        [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds,NomTypeNc]=...
177                            find_pair_indices(PairCiv2,i1_series{2},j1_series{2},MinIndex_i(2),MaxIndex_i(2),MinIndex_j(2),MaxIndex_j(2));
178                    end
179                end
180            case 'displacement'
181                i1_series_Civ1=Param.ActionInput.OriginIndex*ones(size(i1_series{1}));
182                i2_series_Civ1=i1_series{1};i2_series_Civ2=i1_series{1};
183                j1_series_Civ1=ones(size(i1_series{1}));% first j index is 1
184                if isempty(j1_series_Civ1)
185                    j2_series_Civ1=ones(size(i1_series{1}));
186                else
187                    j2_series_Civ1=j1_series_Civ1;
188                end
189                i1_series_Civ2=i1_series_Civ1;
190                j1_series_Civ2=j1_series_Civ1;
191                j2_series_Civ2=j2_series_Civ1;
192                NomTypeNc=Param.InputTable{1,4};
193            case 'PIV volume'
194                % TODO, TODO
195        end
196        if isempty(j1_series_Civ1)
197            FrameIndex_A_Civ1=i1_series_Civ1;
198            FrameIndex_B_Civ1=i2_series_Civ1;
199            j1_series_Civ1=ones(size(i1_series_Civ1));
200            j2_series_Civ1=ones(size(i1_series_Civ1));
201        else
202            FrameIndex_A_Civ1=j1_series_Civ1;
203            FrameIndex_B_Civ1=j2_series_Civ1;
204        end
205        if isempty(PairCiv2)
206            FrameIndex_A_Civ2=FrameIndex_A_Civ1;
207            FrameIndex_B_Civ2=FrameIndex_B_Civ1;
208        else
209            if isempty(j1_series_Civ2)
210                FrameIndex_A_Civ2=i1_series_Civ2;
211                FrameIndex_B_Civ2=i2_series_Civ2;
212                j1_series_Civ2=ones(size(i1_series_Civ2));
213                j2_series_Civ2=ones(size(i1_series_Civ2));
214            else
215                FrameIndex_A_Civ2=j1_series_Civ2;
216                FrameIndex_B_Civ2=j2_series_Civ2;
217            end
218        end
219        if isempty(i1_series_Civ1)||(~isempty(PairCiv2) && isempty(i1_series_Civ2))
220            disp_uvmat('ERROR','no image pair for civ in the input file index range',checkrun)
221            return
222        end
223    end
224   
225    %% check the first image pair
226    try
227        if Param.ActionInput.CheckCiv1% Civ1 is performed
228            first_ima=1;
229            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(1),[],j1_series_Civ1(1));
230            if checkrun
231                if ~exist(ImageName_A,'file')
232                disp_uvmat('ERROR',['first input image ' ImageName_A ' does not exist'],checkrun)
233                return
234                end
235            else
236%                 for ifile=1:numel(i1_series_Civ1)
237%                     if
238            end
239            [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A);
240            FileType_A=FileInfo_A.FileType;
241            if strcmp(FileInfo_A.FileType,'netcdf')
242                FieldName_A=Param.InputFields.FieldName;
243                [DataIn,tild,tild,errormsg]=nc2struct(ImageName_A,{FieldName_A});
244                par_civ1.ImageA=DataIn.(FieldName_A);
245            else
246                [par_civ1.ImageA,VideoObject_A] = read_image(ImageName_A,FileType_A,VideoObject_A,FrameIndex_A_Civ1(1));
247            end
248            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(1),[],j2_series_Civ1(1));
249            if ~exist(ImageName_B,'file')
250                disp_uvmat('ERROR',['first input image ' ImageName_B ' does not exist'],checkrun)
251                return
252            end
253            [FileInfo_B,VideoObject_B]=get_file_info(ImageName_B);
254            FileType_B=FileInfo_B.FileType;
255            if strcmp(FileInfo_B.FileType,'netcdf')
256                FieldName_B=Param.InputFields.FieldName;
257                [DataIn,tild,tild,errormsg]=nc2struct(ImageName_B,{FieldName_B});
258                par_civ1.ImageB=DataIn.(FieldName_B);
259            else
260                [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(1));
261            end
262            NbField=numel(i1_series_Civ1);
263        elseif Param.ActionInput.CheckCiv2 % Civ2 is performed without Civ1
264            ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ2(1),[],j1_series_Civ2(1));
265            if ~exist(ImageName_A,'file')
266                disp_uvmat('ERROR',['first input image ' ImageName_A ' does not exist'],checkrun)
267                return
268            end
269            [FileInfo_A,VideoObject_A]=get_file_info(ImageName_A);
270            FileType_A=FileInfo_A.FileType;
271            [par_civ1.ImageA,VideoObject_A] = read_image(ImageName_A,FileInfo_A.FileType,VideoObject_A,FrameIndex_A_Civ2(1));
272            ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ2(1),[],j2_series_Civ2(1));
273            if ~exist(ImageName_B,'file')
274                disp_uvmat('ERROR',['first input image ' ImageName_B ' does not exist'],checkrun)
275                return
276            end
277            [FileInfo_B,VideoObject_B]=get_file_info(ImageName_B);
278            FileType_B=FileInfo_B.FileType;
279            [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ2(1));
280            NbField=numel(i1_series_Civ2);
281        else
282            NbField=numel(i1_series_Civ1);% no image used (only fix or patch) TO CHECK
283        end
284    catch ME
285        if ~isempty(ME.message)
286            disp_uvmat('ERROR', ['error reading input image: ' ME.message],checkrun)
287            return
288        end
289    end
290   
291   
292    %% Output directory
293    OutputDir='';
294    if CheckOutputFile
295        OutputDir=[Param.OutputSubDir Param.OutputDirExt];
296    end
297end
298
299%% prepare output Data
300ListGlobalAttribute={'Conventions','Program','CivStage'};
301Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
302Data.Program='civ_series';
303Data.CivStage=0;%default
304check_civx=0;%default
305
306%% get timing from the ImaDoc file or input video
307if iview_A~=0
308    XmlFileName=find_imadoc(RootPath_A,SubDir_A,RootFile_A,FileExt_A);
309    Time=[];
310    if ~isempty(XmlFileName)
311        XmlData=imadoc2struct(XmlFileName);
312        if isfield(XmlData,'Time')
313            Time=XmlData.Time;
314            TimeSource='xml';
315        end
316        if isfield(XmlData,'Camera')
317            if isfield(XmlData.Camera,'NbSlice')&& ~isempty(XmlData.Camera.NbSlice)
318                NbSlice_calib{iview}=XmlData.Camera.NbSlice;% Nbre of slices for Zindex in phys transform
319                if ~isequal(NbSlice_calib{iview},NbSlice_calib{1})
320                    msgbox_uvmat('WARNING','inconsistent number of Z indices for the two field series');
321                end
322            end
323            if isfield(XmlData.Camera,'TimeUnit')&& ~isempty(XmlData.Camera.TimeUnit)
324                TimeUnit=XmlData.Camera.TimeUnit;
325            end
326        end
327    end
328    if isempty(Time) && ~isempty(find(strcmp(FileType_A,{'mmreader','video','cine_phantom'})))% case of video input
329        Time=zeros(FileInfo_A.NumberOfFrames+1,2);
330        Time(:,2)=(0:1/FileInfo_A.FrameRate:(FileInfo_A.NumberOfFrames)/FileInfo_A.FrameRate)';
331        TimeSource='video';
332        ColorType='truecolor';
333    end
334    if isempty(Time)% Time = index i +0.001 index j by default
335        %MinIndex_i=min(i1_series_Civ1);
336        MaxIndex_i=max(i2_series_Civ1);
337        %MinIndex_j=min(j1_series_Civ1);
338        MaxIndex_j=max(j2_series_Civ1);
339        Time=(1:MaxIndex_i)'*ones(1,MaxIndex_j);
340        Time=Time+0.001*ones(MaxIndex_i,1)*(1:MaxIndex_j);
341        Time=[zeros(1,MaxIndex_j);Time];% insert a first line of zeros
342        Time=[zeros(MaxIndex_i+1,1) Time];% insert a first column of zeros
343    end
344   
345    if length(FileInfo_A) >1 %case of image with multiple frames
346        nbfield=length(FileInfo_A);
347        nbfield_j=1;
348    end
349end
350
351%%%%% MAIN LOOP %%%%%%
352maskname='';% initiate the mask name
353tic;
354CheckOverwrite=1;%default
355if isfield(Param,'CheckOverwrite')
356    CheckOverwrite=Param.CheckOverwrite;
357end
358for ifield=1:NbField
359    if ~isempty(RUNHandle)% update the waitbar in interactive mode with GUI series  (checkrun=1)
360        update_waitbar(WaitbarHandle,ifield/NbField)
361        if  checkrun && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
362            disp('program stopped by user')
363            break
364        end
365    end
366    if CheckInputFile
367        if iview_A==0 % no nc file has been entered
368            ncfile=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},Param.InputTable{1,5},...
369                NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),j1_series_Civ1(ifield),j2_series_Civ1(ifield));
370        else% an existing nc file has been entered
371            if iview_A==1% if Civ1 is performed
372                Civ1Dir=OutputDir;
373            else
374                Civ1Dir=Param.InputTable{1,2};
375            end
376            if strcmp(Param.ActionInput.ListCompareMode,'PIV')
377                ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),...
378                    j1_series_Civ1(ifield),j2_series_Civ1(ifield));
379            else
380                ncfile=fullfile_uvmat(RootPath_A,Civ1Dir,RootFile_A,'.nc',NomTypeNc,i2_series_Civ1(ifield),[],...
381                    j1_series_Civ1(ifield),j2_series_Civ1(ifield));
382            end
383        end
384        ncfile_out=ncfile;% by default
385        if isfield (Param.ActionInput,'Civ2')
386            i1_civ2=i1_series_Civ2(ifield);
387            i2_civ2=i1_civ2;
388            if ~isempty(i2_series_Civ2)
389                i2_civ2=i2_series_Civ2(ifield);
390            end
391            j1_civ2=1;
392            if ~isempty(j1_series_Civ2)
393                j1_civ2=j1_series_Civ2(ifield);
394            end
395            j2_civ2=i1_civ2;
396            if ~isempty(j2_series_Civ2)
397                j2_civ2=j2_series_Civ2(ifield);
398            end
399            if strcmp(Param.ActionInput.ListCompareMode,'PIV')
400                ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1_civ2,i2_civ2,j1_civ2,j2_civ2);
401            else % displacement
402                ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i2_civ2,[],j2_civ2);
403            end
404        end
405        if ~CheckOverwrite && exist(ncfile_out,'file')
406            disp(['existing output file ' ncfile_out ' already exists, skip to next field'])
407            continue% skip iteration if the mode overwrite is desactivated and the result file already exists
408        end
409    end
410    %% Civ1
411    % if Civ1 computation is requested
412    if isfield (Param.ActionInput,'Civ1')
413        if CheckInputFile
414        disp('civ1 started')
415        end
416        par_civ1=Param.ActionInput.Civ1;
417        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
418            try
419                ImageName_A=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
420                if strcmp(FileInfo_A.FileType,'netcdf')% case of input images in format netcdf
421                    FieldName_A=Param.InputFields.FieldName;
422                    [DataIn,tild,tild,errormsg]=nc2struct(ImageName_A,{FieldName_A});
423                    par_civ1.ImageA=DataIn.(FieldName_A);
424                else % usual image formats for image A
425                    [par_civ1.ImageA,VideoObject_A] = read_image(ImageName_A,FileType_A,VideoObject_A,FrameIndex_A_Civ1(ifield));
426                end
427                ImageName_B=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield));
428                if strcmp(FileInfo_B.FileType,'netcdf') % case of input images in format netcdf
429                    FieldName_B=Param.InputFields.FieldName;
430                    [DataIn,tild,tild,errormsg]=nc2struct(ImageName_B,{FieldName_B});
431                    par_civ1.ImageB=DataIn.(FieldName_B);
432                else % usual image formats for image B
433                    [par_civ1.ImageB,VideoObject_B] = read_image(ImageName_B,FileType_B,VideoObject_B,FrameIndex_B_Civ1(ifield));
434                end
435            catch ME % display errors in reading input images
436                if ~isempty(ME.message)
437                    disp_uvmat('ERROR', ['error reading input image: ' ME.message],checkrun)
438                    continue
439                end
440            end
441            par_civ1.ImageWidth=size(par_civ1.ImageA,2);%FileInfo_A.Width;
442            par_civ1.ImageHeight=size(par_civ1.ImageA,1);%FileInfo_A.Height;
443            list_param=(fieldnames(Param.ActionInput.Civ1))';
444            list_param(strcmp('TestCiv1',list_param))=[];% remove the parameter TestCiv1 from the list
445            Civ1_param=regexprep(list_param,'^.+','Civ1_$0');% insert 'Civ1_' before  each string in list_param
446            Civ1_param=[{'Civ1_ImageA','Civ1_ImageB','Civ1_Time','Civ1_Dt'} Civ1_param]; %insert the names of the two input images
447            %indicate the values of all the global attributes in the output data
448            Data.Civ1_ImageA=ImageName_A;
449            Data.Civ1_ImageB=ImageName_B;
450            i1=i1_series_Civ1(ifield);
451            i2=i1;
452            if ~isempty(i2_series_Civ1)
453                i2=i2_series_Civ1(ifield);
454            end
455            j1=1;
456            if ~isempty(j1_series_Civ1)
457                j1=j1_series_Civ1(ifield);
458            end
459            j2=j1;
460            if ~isempty(j2_series_Civ1)
461                j2=j2_series_Civ1(ifield);
462            end
463            if strcmp(Param.ActionInput.ListCompareMode,'displacement')
464                Data.Civ1_Time=Time(i2+1,j2+1);% the Time is the Time of the secodn image
465                Data.Civ1_Dt=1;% Time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
466            else
467            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
468            Data.Civ1_Dt=Time(i2+1,j2+1)-Time(i1+1,j1+1);
469            end
470            for ilist=1:length(list_param)
471                Data.(Civ1_param{4+ilist})=Param.ActionInput.Civ1.(list_param{ilist});
472            end
473            Data.ListGlobalAttribute=[ListGlobalAttribute Civ1_param];
474            Data.CivStage=1;
475        end
476        % set the list of variables
477        Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_F','Civ1_C'};%  cell array containing the names of the fields to record
478        Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
479        Data.VarAttribute{1}.Role='coord_x';
480        Data.VarAttribute{2}.Role='coord_y';
481        Data.VarAttribute{3}.Role='vector_x';
482        Data.VarAttribute{4}.Role='vector_y';
483        Data.VarAttribute{5}.Role='warnflag';
484        if par_civ1.CheckMask&&~isempty(par_civ1.Mask)
485            if strcmp(maskname,par_civ1.Mask)% mask exist, not already read in civ1
486                par_civ1.Mask=mask; %use mask already opened
487            else
488                try
489                    par_civ1.Mask=imread(par_civ1.Mask);%update the mask, an store it for future use
490                catch ME
491                    if ~isempty(ME.message)
492                        errormsg=['error reading input image: ' ME.message];
493                        disp_uvmat('ERROR',errormsg,checkrun)
494                        return
495                    end
496                end
497                mask=par_civ1.Mask;
498                maskname=par_civ1.Mask;
499            end
500        end
501        if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume')
502            Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
503            Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
504            Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];Data.Civ1_F=[];
505            for ivol=1:NbSlice
506                % caluclate velocity data (y and v in indices, reverse to y component)
507                [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
508                if ~isempty(errormsg)
509                    disp_uvmat('ERROR',errormsg,checkrun)
510                    return
511                end
512                Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
513                Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
514                Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates
515                Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)];
516                Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)];
517                Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
518                Data.Civ1_F=[Data.Civ1_C reshape(F,[],1)];
519            end
520        else %usual PIV
521            % caluclate velocity data (y and v in indices, reverse to y component)
522            [xtable, ytable, utable, vtable, ctable, F, result_conv, errormsg] = civ (par_civ1);
523            if ~isempty(errormsg)
524                disp_uvmat('ERROR',errormsg,checkrun)
525                return
526            end
527            Data.Civ1_X=reshape(xtable,[],1);
528            Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1);
529            Data.Civ1_U=reshape(utable,[],1);
530            Data.Civ1_V=reshape(-vtable,[],1);
531            Data.Civ1_C=reshape(ctable,[],1);
532            Data.Civ1_F=reshape(F,[],1);
533        end
534    else% we use existing Civ1 data
535        if exist('ncfile','var')
536            CivFile=ncfile;
537            [Data,tild,tild,errormsg]=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
538            if ~isempty(errormsg)
539                disp_uvmat('ERROR',errormsg,checkrun)
540                return
541            end
542            [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
543        elseif isfield(Param,'Civ1_X')
544            Data.ListGlobalAttribute={};
545            Data.ListVarName={};
546            Data.VarDimName={};
547            Data.Civ1_X=Param.Civ1_X;
548            Data.Civ1_Y=Param.Civ1_Y;
549            Data.Civ1_U=Param.Civ1_U;
550            Data.Civ1_V=Param.Civ1_V;
551            Data.Civ1_FF=Param.Civ1_FF;
552        end
553    end
554   
555    %% Fix1
556    if isfield (Param.ActionInput,'Fix1')
557         disp('fix1 started')
558        if ~isfield (Param.ActionInput,'Civ1')% if we use existing Civ1, remove previous data beyond Civ1
559            Fix1_attr=find(strcmp('Fix1',Data.ListGlobalAttribute));
560            Data.ListGlobalAttribute(Fix1_attr)=[];
561            for ilist=1:numel(Fix1_attr)
562                Data=rmfield(Data,Data.ListGlobalAttribute{Fix1_attr(ilist)});
563            end
564        end
565        list_param=fieldnames(Param.ActionInput.Fix1)';
566        Fix1_param=regexprep(list_param,'^.+','Fix1_$0');% insert 'Fix1_' before  each string in ListFixParam
567        %indicate the values of all the global attributes in the output data
568        for ilist=1:length(list_param)
569            Data.(Fix1_param{ilist})=Param.ActionInput.Fix1.(list_param{ilist});
570        end
571        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix1_param];
572        Data.ListVarName=[Data.ListVarName {'Civ1_FF'}];
573        Data.VarDimName=[Data.VarDimName {'nb_vec_1'}];
574        nbvar=length(Data.ListVarName);
575        Data.VarAttribute{nbvar}.Role='errorflag';
576        Data.Civ1_FF=int8(fix(Param.ActionInput.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V));
577        Data.CivStage=2;
578    end
579    %% Patch1
580    if isfield (Param.ActionInput,'Patch1')
581        disp('patch1 started')
582        if check_civx
583            errormsg='Civ Matlab input needed for patch';
584            disp_uvmat('ERROR',errormsg,checkrun)
585            return
586        end
587       
588        % record the processing parameters of Patch1 as global attributes in the result nc file
589        list_param=fieldnames(Param.ActionInput.Patch1)';
590        list_param(strcmp('TestPatch1',list_param))=[];% remove 'TestPatch1' from the list of parameters
591        Patch1_param=regexprep(list_param,'^.+','Patch1_$0');% insert 'Patch1_' before  each parameter name
592        for ilist=1:length(list_param)
593            Data.(Patch1_param{ilist})=Param.ActionInput.Patch1.(list_param{ilist});
594        end
595        Data.CivStage=3;% record the new state of processing
596        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch1_param];
597       
598        % list the variables to record
599        nbvar=length(Data.ListVarName);
600        Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
601        Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
602            {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
603        Data.VarAttribute{nbvar+1}.Role='vector_x';
604        Data.VarAttribute{nbvar+2}.Role='vector_y';
605        Data.VarAttribute{nbvar+5}.Role='coord_tps';
606        Data.VarAttribute{nbvar+6}.Role='vector_x';
607        Data.VarAttribute{nbvar+7}.Role='vector_y';
608        Data.Civ1_U_smooth=Data.Civ1_U; % zeros(size(Data.Civ1_X));
609        Data.Civ1_V_smooth=Data.Civ1_V; %zeros(size(Data.Civ1_X));
610        if isfield(Data,'Civ1_FF')
611            ind_good=find(Data.Civ1_FF==0);
612        else
613            ind_good=1:numel(Data.Civ1_X);
614        end
615       
616        % perform Patch calculation using the UVMAT fct 'filter_tps'
617        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
618            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);
619        Data.Civ1_U_smooth(ind_good)=Ures;% take the interpolated (smoothed) velocity values for good vectors, keep civ1 data for the other
620        Data.Civ1_V_smooth(ind_good)=Vres;
621        Data.Civ1_FF(ind_good)=int8(FFres);
622        disp('patch1 performed')
623    end
624   
625    %% Civ2
626    if isfield (Param.ActionInput,'Civ2')
627        disp('civ2 started')
628        par_civ2=Param.ActionInput.Civ2;
629        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
630            par_civ2.ImageA=[];
631            par_civ2.ImageB=[];
632%             i1_civ2=i1_series_Civ2(ifield);
633%             i2_civ2=i1_civ2;
634%             if ~isempty(i2_series_Civ2)
635%                 i2_civ2=i2_series_Civ2(ifield);
636%             end
637%             j1_civ2=1;
638%             if ~isempty(j1_series_Civ2)
639%                 j1_civ2=j1_series_Civ2(ifield);
640%             end
641%             j2_civ2=i1_civ2;
642%             if ~isempty(j2_series_Civ2)
643%                 j2_civ2=j2_series_Civ2(ifield);
644%             end
645            ImageName_A_Civ2=fullfile_uvmat(RootPath_A,SubDir_A,RootFile_A,FileExt_A,NomType_A,i1_civ2,[],j1_civ2);
646           
647            if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1(ifield),FrameIndex_A_Civ2(ifield))
648                par_civ2.ImageA=par_civ1.ImageA;
649            else
650                [par_civ2.ImageA,VideoObject_A] = read_image(ImageName_A_Civ2,FileType_A,VideoObject_A,FrameIndex_A_Civ2(ifield));
651            end
652            ImageName_B_Civ2=fullfile_uvmat(RootPath_B,SubDir_B,RootFile_B,FileExt_B,NomType_B,i2_civ2,[],j2_civ2);
653            if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1(ifield),FrameIndex_B_Civ2)
654                par_civ2.ImageB=par_civ1.ImageB;
655            else
656                [par_civ2.ImageB,VideoObject_B] = read_image(ImageName_B_Civ2,FileType_B,VideoObject_B,FrameIndex_B_Civ2(ifield));
657            end
658%             if strcmp(Param.ActionInput.ListCompareMode,'PIV')
659%                 ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i1,i2,j1,j2);
660%             else % displacement
661%                 ncfile_out=fullfile_uvmat(RootPath_A,OutputDir,RootFile_A,'.nc',NomTypeNc,i2,[],j2);
662%             end
663            par_civ2.ImageWidth=FileInfo_A.Width;
664            par_civ2.ImageHeight=FileInfo_A.Height;
665            if isfield(par_civ2,'Grid')% grid points set as input file
666                if ischar(par_civ2.Grid)%read the grid file if the input is a file name
667                    par_civ2.Grid=dlmread(par_civ2.Grid);
668                    par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
669                end
670            else% automatic grid
671                minix=floor(par_civ2.Dx/2)-0.5;
672                maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
673                miniy=floor(par_civ2.Dy/2)-0.5;
674                maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
675                [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
676                par_civ2.Grid(:,1)=reshape(GridX,[],1);
677                par_civ2.Grid(:,2)=reshape(GridY,[],1);
678            end
679        end
680       
681        % get the guess from patch1 or patch2 (case 'CheckCiv3')
682        if CheckInputFile % read input images (except in mode Test where it is introduced directly in Param.ActionInput.Civ1.ImageNameA and B)
683            if isfield (par_civ2,'CheckCiv3') && par_civ2.CheckCiv3 %get the guess from  patch2
684                SubRange= Data.Civ2_SubRange;
685                NbCentres=Data.Civ2_NbCentres;
686                Coord_tps=Data.Civ2_Coord_tps;
687                U_tps=Data.Civ2_U_tps;
688                V_tps=Data.Civ2_V_tps;
689                CivStage=Data.CivStage;%store the current CivStage
690                Civ1_Dt=Data.Civ2_Dt;
691                Data=[];%reinitialise the result structure Data
692                Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
693                Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
694                Data.Program='civ_series';
695                Data.CivStage=CivStage+1;%update the current civStage after reinitialisation of Data
696                Data.ListVarName={};
697                Data.VarDimName={};
698            else % get the guess from patch1
699                SubRange= Data.Civ1_SubRange;
700                NbCentres=Data.Civ1_NbCentres;
701                Coord_tps=Data.Civ1_Coord_tps;
702                U_tps=Data.Civ1_U_tps;
703                V_tps=Data.Civ1_V_tps;
704                Civ1_Dt=Data.Civ1_Dt;
705                Data.CivStage=4;
706            end
707        else
708            SubRange= par_civ2.Civ1_SubRange;
709            NbCentres=par_civ2.Civ1_NbCentres;
710            Coord_tps=par_civ2.Civ1_Coord_tps;
711            U_tps=par_civ2.Civ1_U_tps;
712            V_tps=par_civ2.Civ1_V_tps;
713            Civ1_Dt=par_civ2.Civ1_Dt;
714            Civ2_Dt=par_civ2.Civ1_Dt;
715            Data.ListVarName={};
716            Data.VarDimName={};
717        end
718        Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
719        Shifty=zeros(size(par_civ2.Grid,1),1);
720        nbval=zeros(size(par_civ2.Grid,1),1);% nbre of interpolated values at each grid point (from the different patch subdomains)
721        if par_civ2.CheckDeformation
722            DUDX=zeros(size(par_civ2.Grid,1),1);
723            DUDY=zeros(size(par_civ2.Grid,1),1);
724            DVDX=zeros(size(par_civ2.Grid,1),1);
725            DVDY=zeros(size(par_civ2.Grid,1),1);
726        end
727        NbSubDomain=size(SubRange,3);
728        for isub=1:NbSubDomain% for each sub-domain of Patch1
729            nbvec_sub=NbCentres(isub);% nbre of Civ vectors in the subdomain
730            ind_sel=find(par_civ2.Grid(:,1)>=SubRange(1,1,isub) & par_civ2.Grid(:,1)<=SubRange(1,2,isub) &...
731                par_civ2.Grid(:,2)>=SubRange(2,1,isub) & par_civ2.Grid(:,2)<=SubRange(2,2,isub));% grid points in the subdomain
732            if ~isempty(ind_sel)
733                epoints = par_civ2.Grid(ind_sel,:);% coordinates of interpolation sites (measurement grids)
734                ctrs=Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
735                nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for each interpolation point (in case of subdomain overlap)
736                EM = tps_eval(epoints,ctrs);% thin plate spline (tps) coefficient
737                Shiftx(ind_sel)=Shiftx(ind_sel)+EM*U_tps(1:nbvec_sub+3,isub);%velocity shift estimated by tps from civ1
738                Shifty(ind_sel)=Shifty(ind_sel)+EM*V_tps(1:nbvec_sub+3,isub);
739                if par_civ2.CheckDeformation
740                    [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
741                    DUDX(ind_sel)=DUDX(ind_sel)+EMDX*U_tps(1:nbvec_sub+3,isub);
742                    DUDY(ind_sel)=DUDY(ind_sel)+EMDY*U_tps(1:nbvec_sub+3,isub);
743                    DVDX(ind_sel)=DVDX(ind_sel)+EMDX*V_tps(1:nbvec_sub+3,isub);
744                    DVDY(ind_sel)=DVDY(ind_sel)+EMDY*V_tps(1:nbvec_sub+3,isub);
745                end
746            end
747        end
748        if par_civ2.CheckMask&&~isempty(par_civ2.Mask)
749            if strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1
750                par_civ2.Mask=mask; %use mask already opened
751            else
752                par_civ2.Mask=imread(par_civ2.Mask);%update the mask, and store it for future use
753                mask=par_civ2.Mask;
754                maskname=par_civ2.Mask;
755            end
756        end
757        if CheckInputFile % else Dt given by par_civ2
758            if strcmp(Param.ActionInput.ListCompareMode,'displacement')
759                Civ1_Dt=1;
760                Civ2_Dt=1;
761            else
762                Civ2_Dt=Time(i2_civ2+1,j2_civ2+1)-Time(i1_civ2+1,j1_civ2+1);
763            end
764        end
765        par_civ2.SearchBoxShift=(Civ2_Dt/Civ1_Dt)*[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
766        % shift the grid points by half the expected shift to provide the correlation box position in image A
767        par_civ2.Grid=[par_civ2.Grid(nbval>=1,1)-par_civ2.SearchBoxShift(:,1)/2 par_civ2.Grid(nbval>=1,2)-par_civ2.SearchBoxShift(:,2)/2];
768        if par_civ2.CheckDeformation
769            par_civ2.DUDX=DUDX(nbval>=1)./nbval(nbval>=1);
770            par_civ2.DUDY=DUDY(nbval>=1)./nbval(nbval>=1);
771            par_civ2.DVDX=DVDX(nbval>=1)./nbval(nbval>=1);
772            par_civ2.DVDY=DVDY(nbval>=1)./nbval(nbval>=1);
773        end
774       
775        % calculate velocity data (y and v in image indices, reverse to y component)
776        [xtable, ytable, utable, vtable, ctable, F,result_conv,errormsg] = civ (par_civ2);
777       
778        list_param=(fieldnames(Param.ActionInput.Civ2))';
779        list_param(strcmp('TestCiv2',list_param))=[];% remove the parameter TestCiv2 from the list
780        Civ2_param=regexprep(list_param,'^.+','Civ2_$0');% insert 'Civ2_' before  each string in list_param
781        Civ2_param=[{'Civ2_ImageA','Civ2_ImageB','Civ2_Time','Civ2_Dt'} Civ2_param]; %insert the names of the two input images
782        %indicate the values of all the global attributes in the output data
783        if exist('ImageName_A','var')
784            Data.Civ2_ImageA=ImageName_A;
785            Data.Civ2_ImageB=ImageName_B;
786             if strcmp(Param.ActionInput.ListCompareMode,'displacement')
787                Data.Civ2_Time=Time(i2_civ2+1,j2_civ2+1);% the Time is the Time of the secodn image
788                Data.Civ2_Dt=1;% Time interval is 1, to yield displacement instead of velocity=displacement/Dt at reading
789             else
790            Data.Civ2_Time=(Time(i2_civ2+1,j2_civ2+1)+Time(i1_civ2+1,j1_civ2+1))/2;
791            Data.Civ2_Dt=Civ2_Dt;
792             end
793        end
794        for ilist=1:length(list_param)
795            Data.(Civ2_param{4+ilist})=Param.ActionInput.Civ2.(list_param{ilist});
796        end
797        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param];
798       
799        nbvar=numel(Data.ListVarName);
800        % define the Civ2 variable (if Civ2 data are not replaced from previous calculation)
801        if isempty(find(strcmp('Civ2_X',Data.ListVarName),1))
802            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
803            Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
804            Data.VarAttribute{nbvar+1}.Role='coord_x';
805            Data.VarAttribute{nbvar+2}.Role='coord_y';
806            Data.VarAttribute{nbvar+3}.Role='vector_x';
807            Data.VarAttribute{nbvar+4}.Role='vector_y';
808            Data.VarAttribute{nbvar+5}.Role='warnflag';
809        end
810        Data.Civ2_X=reshape(xtable,[],1);
811        Data.Civ2_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
812        Data.Civ2_U=reshape(utable,[],1);
813        Data.Civ2_V=reshape(-vtable,[],1);
814        Data.Civ2_C=reshape(ctable,[],1);
815        Data.Civ2_F=reshape(F,[],1);
816        disp('civ2 performed')
817    elseif ~isfield(Data,'ListVarName') % we start there, using existing Civ2 data
818        if exist('ncfile','var')
819            CivFile=ncfile;
820            [Data,tild,tild,errormsg]=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
821            if ~isempty(errormsg)
822                disp_uvmat('ERROR',errormsg,checkrun)
823                return
824            end         
825        elseif isfield(Param,'Civ2_X')% use Civ2 data as input in Param (test mode)
826            Data.ListGlobalAttribute={};
827            Data.ListVarName={};
828            Data.VarDimName={};
829            Data.Civ2_X=Param.Civ2_X;
830            Data.Civ2_Y=Param.Civ2_Y;
831            Data.Civ2_U=Param.Civ2_U;
832            Data.Civ2_V=Param.Civ2_V;
833            Data.Civ2_FF=Param.Civ2_FF;
834        end
835    end
836   
837    %% Fix2
838    if isfield (Param.ActionInput,'Fix2')
839        disp('fix2 started')
840        list_param=fieldnames(Param.ActionInput.Fix2)';
841        Fix2_param=regexprep(list_param,'^.+','Fix2_$0');% insert 'Fix1_' before  each string in ListFixParam
842        %indicate the values of all the global attributes in the output data
843        for ilist=1:length(list_param)
844            Data.(Fix2_param{ilist})=Param.ActionInput.Fix2.(list_param{ilist});
845        end
846        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Fix2_param];
847        %
848        %         ListFixParam=fieldnames(Param.ActionInput.Fix2);
849        %         for ilist=1:length(ListFixParam)
850        %             ParamName=ListFixParam{ilist};
851        %             ListName=['Fix2_' ParamName];
852        %             eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
853        %             eval(['Data.' ListName '=Param.ActionInput.Fix2.' ParamName ';'])
854        %         end
855        if check_civx
856            if ~isfield(Data,'fix2')
857                Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2'];
858                Data.fix2=1;
859                Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}];
860                Data.VarDimName=[Data.VarDimName {'nb_vectors2'}];
861            end
862            Data.vec_FixFlag=fix(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y);
863        else
864            Data.ListVarName=[Data.ListVarName {'Civ2_FF'}];
865            Data.VarDimName=[Data.VarDimName {'nb_vec_2'}];
866            nbvar=length(Data.ListVarName);
867            Data.VarAttribute{nbvar}.Role='errorflag';
868            Data.Civ2_FF=double(fix(Param.ActionInput.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V));
869            Data.CivStage=Data.CivStage+1;
870        end
871    end
872   
873    %% Patch2
874    if isfield (Param.ActionInput,'Patch2')
875        disp('patch2 started')
876        list_param=fieldnames(Param.ActionInput.Patch2)';
877        list_param(strcmp('TestPatch2',list_param))=[];% remove the parameter TestCiv1 from the list
878        Patch2_param=regexprep(list_param,'^.+','Patch2_$0');% insert 'Fix1_' before  each string in ListFixParam
879        %indicate the values of all the global attributes in the output data
880        for ilist=1:length(list_param)
881            Data.(Patch2_param{ilist})=Param.ActionInput.Patch2.(list_param{ilist});
882        end
883        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Patch2_param];
884       
885        nbvar=length(Data.ListVarName);
886        Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
887        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
888            {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
889       
890        Data.VarAttribute{nbvar+1}.Role='vector_x';
891        Data.VarAttribute{nbvar+2}.Role='vector_y';
892        Data.VarAttribute{nbvar+5}.Role='coord_tps';
893        Data.VarAttribute{nbvar+6}.Role='vector_x';
894        Data.VarAttribute{nbvar+7}.Role='vector_y';
895        Data.Civ2_U_smooth=Data.Civ2_U;
896        Data.Civ2_V_smooth=Data.Civ2_V;
897        if isfield(Data,'Civ2_FF')
898            ind_good=find(Data.Civ2_FF==0);
899        else
900            ind_good=1:numel(Data.Civ2_X);
901        end
902        [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures,Vres,tild,FFres]=...
903            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);
904        Data.Civ2_U_smooth(ind_good)=Ures;
905        Data.Civ2_V_smooth(ind_good)=Vres;
906        Data.Civ2_FF(ind_good)=FFres;
907        Data.CivStage=Data.CivStage+1;
908        disp('patch2 performed')
909    end
910   
911    %% write result in a netcdf file if requested
912    if CheckOutputFile
913        errormsg=struct2nc(ncfile_out,Data);
914        if isempty(errormsg)
915            disp([ncfile_out ' written'])
916            %[success,msg] = fileattrib(ncfile_out ,'+w','g');% done in struct2nc
917        else
918            disp(errormsg)
919        end
920        disp(['ellapsed time ' num2str(toc/60,2) ' minutes'])
921    end
922end
923
924
925% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
926%--------------------------------------------------------------------------
927% function [xtable ytable utable vtable typevector] = civ (image1,image2,ibx,iby step, subpixfinder, mask, roi)
928%
929% OUTPUT:
930% xtable: set of x coordinates
931% ytable: set of y coordiantes
932% utable: set of u displacements (along x)
933% vtable: set of v displacements (along y)
934% ctable: max image correlation for each vector
935% typevector: set of flags, =1 for good, =0 for NaN vectors
936%
937%INPUT:
938% par_civ: structure of input parameters, with fields:
939%  .ImageA: first image for correlation (matrix)
940%  .ImageB: second image for correlation(matrix)
941%  .CorrBoxSize: 1,2 vector giving the size of the correlation box in x and y
942%  .SearchBoxSize:  1,2 vector giving the size of the search box in x and y
943%  .SearchBoxShift: 1,2 vector or 2 column matrix (for civ2) giving the shift of the search box in x and y
944%  .CorrSmooth: =1 or 2 determines the choice of the sub-pixel determination of the correlation max
945%  .ImageWidth: nb of pixels of the image in x
946%  .Dx, Dy: mesh for the PIV calculation
947%  .Grid: grid giving the PIV calculation points (alternative to .Dx .Dy): centres of the correlation boxes in Image A
948%  .Mask: name of a mask file or mask image matrix itself
949%  .MinIma: thresholds for image luminosity
950%  .MaxIma
951%  .CheckDeformation=1 for subpixel interpolation and image deformation (linear transform)
952%  .DUDX: matrix of deformation obtained from patch at each grid point
953%  .DUDY
954%  .DVDX:
955%  .DVDY
956
957function [xtable,ytable,utable,vtable,ctable,F,result_conv,errormsg] = civ (par_civ)
958
959%% prepare measurement grid
960if isfield(par_civ,'Grid')% grid points set as input, central positions of the sub-images in image A
961    if ischar(par_civ.Grid)%read the grid file if the input is a file name (grid in x, y image coordinates)
962        par_civ.Grid=dlmread(par_civ.Grid);
963        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
964    end
965    % else par_civ.Grid is already an array, no action here
966else% automatic grid in x, y image coordinates
967    minix=floor(par_civ.Dx/2)-0.5;
968    maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
969    miniy=floor(par_civ.Dy/2)-0.5;% first automatic grid point at half the mesh Dy
970    maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
971    [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
972    par_civ.Grid(:,1)=reshape(GridX,[],1);
973    par_civ.Grid(:,2)=reshape(GridY,[],1);% increases with array index
974end
975nbvec=size(par_civ.Grid,1);
976
977%% prepare correlation and search boxes
978ibx2=floor(par_civ.CorrBoxSize(1)/2);
979iby2=floor(par_civ.CorrBoxSize(2)/2);
980isx2=floor(par_civ.SearchBoxSize(1)/2);
981isy2=floor(par_civ.SearchBoxSize(2)/2);
982shiftx=round(par_civ.SearchBoxShift(:,1));%use the input shift estimate, rounded to the next integer value
983shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
984if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
985    shiftx=shiftx*ones(nbvec,1);
986    shifty=shifty*ones(nbvec,1);
987end
988
989%% Array initialisation and default output  if par_civ.CorrSmooth=0 (just the grid calculated, no civ computation)
990xtable=round(par_civ.Grid(:,1)+0.5)-0.5;
991ytable=round(par_civ.ImageHeight-par_civ.Grid(:,2)+0.5)-0.5;% y index corresponding to the position in image coordiantes
992utable=shiftx;%zeros(nbvec,1);
993vtable=shifty;%zeros(nbvec,1);
994ctable=zeros(nbvec,1);
995F=zeros(nbvec,1);
996result_conv=[];
997errormsg='';
998
999%% prepare mask
1000if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
1001    if strcmp(par_civ.Mask,'all')
1002        return    % get the grid only, no civ calculation
1003    elseif ischar(par_civ.Mask)
1004        par_civ.Mask=imread(par_civ.Mask);% read the mask if not allready done
1005    end
1006end
1007check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
1008check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
1009
1010par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
1011par_civ.ImageB=sum(double(par_civ.ImageB),3);
1012[npy_ima npx_ima]=size(par_civ.ImageA);
1013if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima])
1014    errormsg='image pair with unequal size';
1015    return
1016end
1017
1018%% Apply mask
1019% Convention for mask IDEAS TO IMPLEMENT ?
1020% mask >200 : velocity calculated
1021%  200 >=mask>150;velocity not calculated, interpolation allowed (bad spots)
1022% 150>=mask >100: velocity not calculated, nor interpolated
1023%  100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries)
1024%  20>=mask: velocity=0
1025checkmask=0;
1026MinA=min(min(par_civ.ImageA));
1027%MinB=min(min(par_civ.ImageB));
1028%check_undefined=false(size(par_civ.ImageA));
1029if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
1030    checkmask=1;
1031    if ~isequal(size(par_civ.Mask),[npy_ima npx_ima])
1032        errormsg='mask must be an image with the same size as the images';
1033        return
1034    end
1035    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
1036%     par_civ.ImageA(check_undefined)=0;% put image A to zero (i.e. the min image value) in the undefined  area
1037%     par_civ.ImageB(check_undefined)=0;% put image B to zero (i.e. the min image value) in the undefined  area
1038end
1039
1040%% compute image correlations: MAINLOOP on velocity vectors
1041corrmax=0;
1042sum_square=1;% default
1043mesh=1;% default
1044CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
1045if CheckDeformation
1046    mesh=0.25;%mesh in pixels for subpixel image interpolation (x 4 in each direction)
1047    par_civ.CorrSmooth=2;% use SUBPIX2DGAUSS (take into account more points near the max)
1048end
1049
1050if par_civ.CorrSmooth~=0 % par_civ.CorrSmooth=0 implies no civ computation (just input image and grid points given)
1051    for ivec=1:nbvec
1052        iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
1053        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
1054        F(ivec)=0;
1055        subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
1056        subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
1057        subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
1058        subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
1059        image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
1060        image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
1061        check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
1062        check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
1063        check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
1064        check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
1065        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
1066        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
1067        if checkmask
1068            mask1_crop=ones(numel(subrange1_y),numel(subrange1_x));% default value=1 for mask
1069            mask2_crop=ones(numel(subrange2_y),numel(subrange2_x));% default value=1 for mask
1070            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
1071            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
1072            sizemask=sum(sum(mask1_crop))/(numel(subrange1_y)*numel(subrange1_x));%size of the masked part relative to the correlation sub-image
1073            if sizemask > 1/2% eliminate point if more than half of the correlation box is masked
1074                F(ivec)=3; %
1075                utable(ivec)=0;
1076                vtable(ivec)=0;
1077            else
1078                image1_crop=image1_crop.*~mask1_crop;% put to zero the masked pixels (mask1_crop='true'=1)
1079                image2_crop=image2_crop.*~mask2_crop;
1080                image1_mean=mean(mean(image1_crop))/(1-sizemask);
1081                image2_mean=mean(mean(image2_crop))/(1-sizemask);
1082            end
1083        else
1084            image1_mean=mean(mean(image1_crop));
1085            image2_mean=mean(mean(image2_crop));
1086        end
1087        %threshold on image minimum
1088        if F(ivec)~=3
1089            if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
1090                F(ivec)=3;
1091                %threshold on image maximum
1092            elseif check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
1093                F(ivec)=3;
1094            end
1095            if F(ivec)==3
1096                utable(ivec)=0;
1097                vtable(ivec)=0;
1098            else
1099                %mask
1100                if checkmask
1101                    image1_crop=(image1_crop-image1_mean).*~mask1_crop;%substract the mean, put to zero the masked parts
1102                    image2_crop=(image2_crop-image2_mean).*~mask2_crop;
1103                else
1104                    image1_crop=(image1_crop-image1_mean);
1105                    image2_crop=(image2_crop-image2_mean);
1106                end
1107                %deformation
1108                if CheckDeformation
1109                    xi=(1:mesh:size(image1_crop,2));
1110                    yi=(1:mesh:size(image1_crop,1))';
1111                    [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
1112                    XIant=XI-par_civ.DUDX(ivec)*XI+par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
1113                    YIant=YI+par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
1114                    image1_crop=interp2(image1_crop,XIant,YIant);
1115                    image1_crop(isnan(image1_crop))=0;
1116                    xi=(1:mesh:size(image2_crop,2));
1117                    yi=(1:mesh:size(image2_crop,1))';
1118                    image2_crop=interp2(image2_crop,xi,yi,'*spline');
1119                    image2_crop(isnan(image2_crop))=0;
1120                end
1121                sum_square=sum(sum(image1_crop.*image1_crop));
1122                %reference: Oliver Pust, PIV: Direct Cross-Correlation
1123                result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
1124                corrmax= max(max(result_conv));
1125                result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
1126                %Find the correlation max, at 255
1127                [y,x] = find(result_conv==255,1);
1128                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
1129                sum_square=sum_square*sum(sum(subimage2_crop.*subimage2_crop));% product of variances of image 1 and 2
1130                sum_square=sqrt(sum_square);% srt of the variance product to normalise correlation
1131                if ~isempty(y) && ~isempty(x)
1132                    try
1133                        if par_civ.CorrSmooth==1
1134                            [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
1135                        elseif par_civ.CorrSmooth==2
1136                            [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
1137                        else
1138                            [vector,F(ivec)] = quadr_fit(result_conv,x,y);
1139                        end
1140                        utable(ivec)=vector(1)*mesh+shiftx(ivec);
1141                        vtable(ivec)=vector(2)*mesh+shifty(ivec);
1142                        xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
1143                        ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
1144                        iref=round(xtable(ivec)+0.5);% nearest image index for the middle of the vector
1145                        jref=round(ytable(ivec)+0.5);
1146                        % eliminate vectors located in the mask
1147                        if  checkmask && (iref<1 || jref<1 ||iref>npx_ima || jref>npy_ima ||( par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100))
1148                            utable(ivec)=0;
1149                            vtable(ivec)=0;
1150                            F(ivec)=3;
1151                        end
1152                        ctable(ivec)=corrmax/sum_square;% correlation value
1153                    catch ME
1154                        F(ivec)=3;
1155                        disp(ME.message)
1156                    end
1157                else
1158                    F(ivec)=3;
1159                end
1160            end
1161        end
1162    end
1163end
1164result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
1165
1166%------------------------------------------------------------------------
1167% --- Find the maximum of the correlation function after interpolation
1168% OUPUT:
1169% vector = optimum displacement vector with subpixel correction
1170% F =flag: =0 OK
1171%           =-2 , warning: max too close to the edge of the search box (1 pixel margin)
1172% INPUT:
1173% x,y: position of the maximum correlation at integer values
1174
1175function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
1176%------------------------------------------------------------------------
1177% vector=[0 0]; %default
1178F=0;
1179[npy,npx]=size(result_conv);
1180result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1  (=0 by discretisation, to avoid divergence in the log)
1181%the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ??? Israel Institute of Technology
1182%http://urapiv.wordpress.com
1183peaky = y;
1184if y < npy && y > 1
1185    f0 = log(result_conv(y,x));
1186    f1 = log(result_conv(y-1,x));
1187    f2 = log(result_conv(y+1,x));
1188    peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
1189else
1190    F=-2; % warning flag for vector truncated by the limited search box
1191end
1192peakx=x;
1193if x < npx-1 && x > 1
1194    f0 = log(result_conv(y,x));
1195    f1 = log(result_conv(y,x-1));
1196    f2 = log(result_conv(y,x+1));
1197    peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
1198else
1199    F=-2; % warning flag for vector truncated by the limited search box
1200end
1201vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
1202
1203%------------------------------------------------------------------------
1204% --- Find the maximum of the correlation function after interpolation
1205function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
1206%------------------------------------------------------------------------
1207% vector=[0 0]; %default
1208F=-2;
1209peaky=y;
1210peakx=x;
1211result_conv(result_conv<1)=1; %set to 1 correlation values smaller than 1 (to avoid divergence in the log)
1212[npy,npx]=size(result_conv);
1213if (x < npx) && (y < npy) && (x > 1) && (y > 1)
1214    F=0;
1215    for i=-1:1
1216        for j=-1:1
1217            %following 15 lines based on
1218            %H. Nobach ??? M. Honkanen (2005)
1219            %Two-dimensional Gaussian regression for sub-pixel displacement
1220            %estimation in particle image velocimetry or particle position
1221            %estimation in particle tracking velocimetry
1222            %Experiments in Fluids (2005) 38: 511???515
1223            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
1224            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
1225            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
1226            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
1227            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
1228        end
1229    end
1230    c10=(1/6)*sum(sum(c10));
1231    c01=(1/6)*sum(sum(c01));
1232    c11=(1/4)*sum(sum(c11));
1233    c20=(1/6)*sum(sum(c20));
1234    c02=(1/6)*sum(sum(c02));
1235    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
1236    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
1237    if abs(deltax)<1
1238        peakx=x+deltax;
1239    end
1240    if abs(deltay)<1
1241        peaky=y+deltay;
1242    end
1243end
1244vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
1245
1246%------------------------------------------------------------------------
1247% --- Find the maximum of the correlation function after quadratic interpolation
1248function [vector,F] = quadr_fit(result_conv,x,y)
1249[npy,npx]=size(result_conv);
1250if x<4 || y<4 || npx-x<4 ||npy-y <4
1251    F=-2;
1252    vector=[x y];
1253else
1254    F=0;
1255    x_ind=x-4:x+4;
1256    y_ind=y-4:y+4;
1257    x_vec=0.25*(x_ind-x);
1258    y_vec=0.25*(y_ind-y);
1259    [X,Y]=meshgrid(x_vec,y_vec);
1260    coord=[reshape(X,[],1) reshape(Y,[],1)];
1261    result_conv=reshape(result_conv(y_ind,x_ind),[],1);
1262   
1263   
1264    % n=numel(X);
1265    % x=[X Y];
1266    % X=X-0.5;
1267    % Y=Y+0.5;
1268    % y = (X.*X+2*Y.*Y+X.*Y+6) + 0.1*rand(n,1);
1269    p = polyfitn(coord,result_conv,2);
1270    A(1,1)=2*p.Coefficients(1);
1271    A(1,2)=p.Coefficients(2);
1272    A(2,1)=p.Coefficients(2);
1273    A(2,2)=2*p.Coefficients(4);
1274    vector=[x y]'-A\[p.Coefficients(3) p.Coefficients(5)]';
1275    vector=vector'-[floor(npx/2) floor(npy/2)]-1 ;
1276    % zg = polyvaln(p,coord);
1277    % figure
1278    % surf(x_vec,y_vec,reshape(zg,9,9))
1279    % hold on
1280    % plot3(X,Y,reshape(result_conv,9,9),'o')
1281    % hold off
1282end
1283
1284%'RUN_FIX': function for fixing velocity fields:
1285%-----------------------------------------------
1286% RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)
1287%
1288%filename: name of the netcdf file (used as input and output)
1289%field: structure specifying the names of the fields to fix (depending on civ1 or civ2)
1290%.vel_type='civ1' or 'civ2';
1291%.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1);
1292%.fixflag=name of fix flag variable ('vec_FixFlag' for civ1)
1293%flagindex: flag specifying which values of vec_f are removed:
1294% if flagindex(1)=1: vec_f=-2 vectors are removed
1295% if flagindex(2)=1: vec_f=3 vectors are removed
1296% if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
1297%iter=1 for civ1 fields and iter=2 for civ2 fields
1298%thresh_vecC: threshold in the image correlation vec_C
1299%flag_mask: =1 mask used to remove vectors (0 else)
1300%maskname: name of the mask image file for fix
1301%thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
1302%inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
1303%fileref: .nc file name for a reference velocity (='': refrence 0 used)
1304%fieldref: 'civ1','filter1'...feld used in fileref
1305
1306function FF=fix(Param,F,C,U,V,X,Y)
1307FF=zeros(size(F));%default
1308
1309%criterium on warn flags
1310FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'};
1311FlagVal=[-2 2 3 4];
1312for iflag=1:numel(FlagName)
1313    if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag})
1314        FF=(FF==1| F==FlagVal(iflag));
1315    end
1316end
1317%criterium on correlation values
1318if isfield (Param,'MinCorr')
1319    FF=FF==1 | C<Param.MinCorr;
1320end
1321if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
1322    Umod= U.*U+V.*V;
1323    if isfield (Param,'MinVel')&&~isempty(Param.MinVel)
1324        FF=FF==1 | Umod<(Param.MinVel*Param.MinVel);
1325    end
1326    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
1327        FF=FF==1 | Umod>(Param.MaxVel*Param.MaxVel);
1328    end
1329end
1330
1331
1332%------------------------------------------------------------------------
1333% --- determine the list of index pairs of processing file
1334function [i1_series,i2_series,j1_series,j2_series,check_bounds,NomTypeNc]=...
1335    find_pair_indices(str_civ,i_series,j_series,MinIndex_i,MaxIndex_i,MinIndex_j,MaxIndex_j)
1336%------------------------------------------------------------------------
1337i1_series=i_series;% set of first image indexes
1338i2_series=i_series;
1339j1_series=j_series;%ones(size(i_series));% set of first image numbers
1340j2_series=j_series;%ones(size(i_series));
1341r=regexp(str_civ,'^\D(?<ind>[i|j])=( -| )(?<num1>\d+)\|(?<num2>\d+)','names');
1342if ~isempty(r)
1343    mode=['D' r.ind];
1344    ind1=str2num(r.num1);
1345    ind2=str2num(r.num2);
1346else
1347    mode='j1-j2';
1348    r=regexp(str_civ,'^j= (?<num1>[a-z])-(?<num2>[a-z])','names');
1349    if ~isempty(r)
1350        NomTypeNc='_1ab';
1351    else
1352        r=regexp(str_civ,'^j= (?<num1>[A-Z])-(?<num2>[A-Z])','names');
1353        if ~isempty(r)
1354            NomTypeNc='_1AB';
1355        else
1356            r=regexp(str_civ,'^j= (?<num1>\d+)-(?<num2>\d+)','names');
1357            if ~isempty(r)
1358                NomTypeNc='_1_1-2';
1359            end
1360        end
1361    end
1362    if isempty(r)
1363        display('wrong pair mode input option')
1364    else
1365        ind1=stra2num(r.num1);
1366        ind2=stra2num(r.num2);
1367    end
1368end
1369switch mode
1370    case 'Di'
1371        i1_series=i_series-ind1;% set of first image numbers
1372        i2_series=i_series+ind2;
1373        check_bounds=i1_series<MinIndex_i | i2_series>MaxIndex_i;
1374        if isempty(j_series)
1375            NomTypeNc='_1-2';
1376        else
1377            j1_series=j_series;
1378            j2_series=j_series;
1379            NomTypeNc='_1-2_1';
1380        end
1381    case 'Dj'
1382        j1_series=j_series-ind1;
1383        j2_series=j_series+ind2;
1384        check_bounds=j1_series<MinIndex_j | j2_series>MaxIndex_j;
1385        NomTypeNc='_1_1-2';
1386    otherwise %bursts
1387        i1_series=i_series(1,:);% do not sweep the j index
1388        i2_series=i_series(1,:);
1389        j1_series=ind1*ones(1,size(i_series,2));% j index is fixed by pair choice
1390        j2_series=ind2*ones(1,size(i_series,2));
1391        check_bounds=zeros(size(i1_series));% no limitations due to min-max indices
1392end
1393
1394
1395
1396
Note: See TracBrowser for help on using the repository browser.