Ignore:
Timestamp:
Apr 2, 2013, 9:13:42 AM (11 years ago)
Author:
sommeria
Message:

various bugs repaired . civ_series further developed

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/series/civ_series.m

    r597 r598  
    1 %'merge_proj': concatene several fields from series, can project them on a regular grid in phys coordinates
     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
    26%------------------------------------------------------------------------
    3 % function ParamOut=merge_proj(Param)
    4 %------------------------------------------------------------------------
    5 %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
     7% function [Data,errormsg,result_conv]= civ_series(Param,ncfile)
    68%
    79%OUTPUT
    8 % ParamOut: sets options in the GUI series.fig needed for the function
     10% Data=structure containing the PIV results and information on the processing parameters
     11% errormsg=error message char string, default=''
     12% resul_conv: image inter-correlation function for the last grid point (used for tests)
    913%
    1014%INPUT:
    11 % In run mode, the input parameters are given as a Matlab structure Param copied from the GUI series.
    12 % In batch mode, Param is the name of the corresponding xml file containing the same information
    13 % when Param.Action.RUN=0 (as activated when the current Action is selected
    14 % in series), the function ouput paramOut set the activation of the needed GUI elements
     15% Param: input images and processing parameters
     16%     .Civ1: for civ1
     17%     .Fix1:
     18%     .Patch1:
     19%     .Civ2: for civ2
     20%     .Fix2:
     21%     .Patch2:
     22% ncfile: name of a netcdf file to be created for the result (extension .nc)
    1523%
    16 % Param contains the elements:(use the menu bar command 'export/GUI config' in series to
    17 % see the current structure Param)
    18 %    .InputTable: cell of input file names, (several lines for multiple input)
    19 %                      each line decomposed as {RootPath,SubDir,Rootfile,NomType,Extension}
    20 %    .OutputSubDir: name of the subdirectory for data outputs
    21 %    .OutputDirExt: directory extension for data outputs
    22 %    .Action: .ActionName: name of the current activated function
    23 %             .ActionPath:   path of the current activated function
    24 %             .ActionExt: fct extension ('.m', Matlab fct, '.sh', compiled   Matlab fct
    25 %             .RUN =0 for GUI input, =1 for function activation
    26 %             .RunMode='local','background', 'cluster': type of function  use
    27 %             
    28 %    .IndexRange: set the file or frame indices on which the action must be performed
    29 %    .FieldTransform: .TransformName: name of the selected transform function
    30 %                     .TransformPath:   path  of the selected transform function
    31 %    .InputFields: sub structure describing the input fields withfields
    32 %              .FieldName: name(s) of the field
    33 %              .VelType: velocity type
    34 %              .FieldName_1: name of the second field in case of two input series
    35 %              .VelType_1: velocity type of the second field in case of two input series
    36 %              .Coord_y: name of y coordinate variable
    37 %              .Coord_x: name of x coordinate variable
    38 %    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
    39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%civ_series
    40 
    41 function ParamOut=civ_series(Param)
    42 
     24%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     25%  Copyright  2011, LEGI / CNRS-UJF-INPG, joel.sommeria@legi.grenoble-inp.fr.
     26%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     27%     This is part of the toolbox UVMAT.
     28%
     29%     UVMAT is free software; you can redistribute it and/or modify
     30%     it under the terms of the GNU General Public License as published by
     31%     the Free Software Foundation; either version 2 of the License, or
     32%     (at your option) any later version.
     33%
     34%     UVMAT is distributed in the hope that it will be useful,
     35%     but WITHOUT ANY WARRANTY; without even the implied warranty of
     36%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     37%     GNU General Public License (open UVMAT/COPYING.txt) for more details.
     38%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
     39
     40function [Data,errormsg,result_conv]= civ_series(Param,ncfile)
     41errormsg='';
     42path_series=fileparts(which('series'));
     43addpath(fullfile(path_series,'series'))
    4344%% set the input elements needed on the GUI series when the action is selected in the menu ActionName
    4445if isstruct(Param) && isequal(Param.Action.RUN,0)
    45     ParamOut.AllowInputSort='off';...% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    46     ParamOut.WholeIndexRange='off';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
    47     ParamOut.NbSlice='off'; ...%nbre of slices ('off' by default)
    48     ParamOut.VelType='off';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    49     ParamOut.FieldName='off';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    50     ParamOut.FieldTransform = 'off';...%can use a transform function
    51     ParamOut.ProjObject='off';...%can use projection object(option 'off'/'on',
    52     ParamOut.Mask='off';...%can use mask option   (option 'off'/'on', 'off' by default)
    53     ParamOut.OutputDirExt='.civ';%set the output dir extension
    54 return
     46    Data=civ_input(Param);% introduce the civ parameters using the GUI civ_input
     47    Data.Program=mfilename;
     48    Data.AllowInputSort='off';...% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
     49        Data.WholeIndexRange='off';...% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
     50        Data.NbSlice='off'; ...%nbre of slices ('off' by default)
     51        Data.VelType='off';...% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
     52        Data.FieldName='off';...% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
     53        Data.FieldTransform = 'off';...%can use a transform function
     54        Data.ProjObject='off';...%can use projection object(option 'off'/'on',
     55        Data.Mask='off';...%can use mask option   (option 'off'/'on', 'off' by default)
     56        Data.OutputDirExt='.civ';%set the output dir extension
     57    return
    5558end
    5659
    … …  
    6366    checkrun=0;
    6467end
    65 
    66 ParamOut=Param; %default output
    67 if ~isfield(Param,'InputFields')
    68     Param.InputFields.FieldName='';
    69 end
    70 OutputSubDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
    71 
    72 %% root input file type
    73 RootPath=Param.InputTable(:,1);
    74 RootFile=Param.InputTable(:,3);
    75 SubDir=Param.InputTable(:,2);
    76 NomType=Param.InputTable(:,4);
    77 FileExt=Param.InputTable(:,5);
    78 [filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
    79 %%%%%%%%%%%%
    80 % The cell array filecell is the list of input file names, while
    81 % filecell{iview,fileindex}:
    82 %        iview: line in the table corresponding to a given file series
    83 %        fileindex: file index within  the file series,
    84 % i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
    85 % i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
    86 %%%%%%%%%%%%
    87 NbSlice=1;%default
    88 if isfield(Param.IndexRange,'NbSlice')&&~isempty(Param.IndexRange.NbSlice)
    89     NbSlice=Param.IndexRange.NbSlice;
    90 end
    91 nbview=numel(i1_series);%number of input file series (lines in InputTable)
    92 nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
    93 nbfield_i=size(i1_series{1},2); %nb of fields for the i index
    94 nbfield=nbfield_j*nbfield_i; %total number of fields
    95 nbfield_i=floor(nbfield/NbSlice);%total number of  indexes in a slice (adjusted to an integer number of slices)
    96 nbfield=nbfield_i*NbSlice; %total number of fields after adjustement
    97 
    98 %determine the file type on each line from the first input file
    99 ImageTypeOptions={'image','multimage','mmreader','video'};
    100 NcTypeOptions={'netcdf','civx','civdata'};
    101 for iview=1:nbview
    102     if ~exist(filecell{iview,1}','file')
    103         displ_uvmat('ERROR',['the first input file ' filecell{iview,1} ' does not exist'],checkrun)
    104         return
    105     end
    106     [FileType{iview},FileInfo{iview},MovieObject{iview}]=get_file_type(filecell{iview,1});
    107     CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images
    108     CheckNc{iview}=~isempty(find(strcmp(FileType{iview},NcTypeOptions)));% =1 for netcdf files
    109     if ~isempty(j1_series{iview})
    110         frame_index{iview}=j1_series{iview};
     68OutputDir=[Param.OutputSubDir Param.OutputDirExt];
     69
     70Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
     71Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
     72Data.Program='civ_series';
     73Data.CivStage=0;%default
     74ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'}; %variables to read
     75ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'};
     76mask='';
     77maskname='';%default
     78check_civx=0;%default
     79check_civ1=0;%default
     80check_patch1=0;%default
     81
     82% case of input Param set by an xml file (batch mode)
     83if ischar(Param)
     84    Param=xml2struct(Param); %if Param is the name of an xml file, read this file as a Matlab structure
     85end
     86
     87RootPath=Param.InputTable{1,1};
     88RootFile=Param.InputTable{1,3};
     89SubDir=Param.InputTable{1,2};
     90NomType=Param.InputTable{1,4};
     91FileExt=Param.InputTable{1,5};
     92PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1;
     93PairCiv2='';
     94if isfield(Param.ActionInput.PairIndices,'ListPairCiv2')
     95    PairCiv2=Param.ActionInput.PairIndices.ListPairCiv2;
     96end
     97
     98% option use with GUI series
     99NbField=1;
     100MovieObject_A=[];
     101if isfield(Param,'InputTable')
     102    MaxIndex=cell2mat(Param.IndexRange.MaxIndex);
     103    MinIndex=cell2mat(Param.IndexRange.MinIndex);
     104    [filecell,i_series,tild,j_series]=get_file_series(Param);
     105    [i1_series_Civ1,i2_series_Civ1,j1_series_Civ1,j2_series_Civ1,check_bounds,NomTypeNc]=...
     106        find_pair_indices(PairCiv1,i_series{1},j_series{1},MinIndex,MaxIndex);
     107    if ~isempty(PairCiv2)
     108        [i1_series_Civ2,i2_series_Civ2,j1_series_Civ2,j2_series_Civ2,check_bounds_Civ2]=...
     109            find_pair_indices(PairCiv2,i_series{1},j_series{1},MinIndex,MaxIndex);
     110        check_bounds=check_bounds | check_bounds_Civ2;
     111    end
     112    i1_series_Civ1=i1_series_Civ1(~check_bounds);
     113    i2_series_Civ1=i2_series_Civ1(~check_bounds);
     114    j1_series_Civ1=j1_series_Civ1(~check_bounds);
     115    j2_series_Civ1=j2_series_Civ1(~check_bounds);
     116    if ~isempty(j1_series_Civ1)
     117        FrameIndex_A_Civ1=j1_series_Civ1;
     118        FrameIndex_B_Civ1=j2_series_Civ1;
    111119    else
    112         frame_index{iview}=i1_series{iview};
    113     end
    114 end
    115 
    116 
    117 %% calibration data and timing: read the ImaDoc files
    118 [XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
    119 if size(time,1)>1
    120     diff_time=max(max(diff(time)));
    121     if diff_time>0
    122         displ_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)],checkrun)
    123     end   
    124 end
    125 
    126 %% coordinate transform or other user defined transform
    127 
    128 %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
    129  % EDIT FROM HERE
    130 
    131 %% check the validity of  input file types
    132 if CheckImage{1}
    133     FileExtOut='.png'; % write result as .png images for image inputs
    134 elseif CheckNc{1}
    135     FileExtOut='.nc';% write result as .nc files for netcdf inputs
    136 else
    137     displ_uvmat('ERROR',['invalid file type input ' FileType{1}],checkrun)
    138     return
    139 end
    140 for iview=1:nbview
    141         if ~isequal(CheckImage{iview},CheckImage{1})||~isequal(CheckNc{iview},CheckNc{1})
    142         displ_uvmat('ERROR','input set of input series: need  either netcdf either image series',checkrun)
    143     return
    144     end
    145 end
    146 NomTypeOut=NomType;% output file index will indicate the first and last ref index in the series
    147 % if checkrun==1
    148 %     ParamOut.Specific=[];%no specific parameter
    149 %     return %stop here for interactive input (option Param.Specific='?')
    150 % end
    151 
    152 %% Set field names and velocity types
    153 %use Param.InputFields for all views
    154 
    155 %% MAIN LOOP ON SLICES
    156 %%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
    157 
    158 index_slice=i_slice:NbSlice:nbfield;% select file indices of the slice
    159 nbfiles=0;
    160 nbmissing=0;
    161 
    162 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    163 %TODO: call civ_matlab
    164 for index=index_slice
    165     if checkrun
    166         stopstate=get(Param.RUNHandle,'BusyAction');
    167         update_waitbar(Param.WaitbarHandle,index/nbfield)
    168     else
    169         stopstate='queue';
    170     end
    171     if ~isequal(stopstate,'queue')% enable STOP command
    172         return
    173     end
    174     %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    175     Data=cell(1,nbview);%initiate the set Data
    176     nbtime=0;
    177     for iview=1:nbview
    178         %% reading input file(s)
    179         [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},Param.InputFields,frame_index{iview}(index));
    180         if ~isempty(errormsg)
    181             errormsg=['merge_proj/read_field/' errormsg];
    182             display(errormsg)
    183             break
    184         end
    185         timeread(iview)=0;
    186         if isfield(Data{iview},'Time')
    187             timeread(iview)=Data{iview}.Time;
    188             nbtime=nbtime+1;
    189         end
    190         if ~isempty(NbSlice_calib)
    191             Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
    192         end
     120        FrameIndex_A_Civ1=i1_series_Civ1;
     121        FrameIndex_B_Civ1=i2_series_Civ1;
     122    end
     123    if ~isempty(PairCiv2)
     124        i1_series_Civ2=i1_series_Civ2(~check_bounds);
     125        i2_series_Civ2=i2_series_Civ2(~check_bounds);
     126        j1_series_Civ2=j1_series_Civ2(~check_bounds);
     127        j2_series_Civ2=j2_series_Civ2(~check_bounds);
     128        if ~isempty(j1_series_Civ2)
     129            FrameIndex_A_Civ2=j1_series_Civ2;
     130            FrameIndex_B_Civ2=j2_series_Civ2;
     131        else
     132            FrameIndex_A_Civ2=i1_series_Civ2;
     133            FrameIndex_B_Civ2=i2_series_Civ2;
     134        end
     135    end
     136   
     137    NbField=numel(i1_series_Civ1);
     138    ImageTypeOptions={'image','multimage','mmreader','video'};
     139    [FileType_A,FileInfo,MovieObject_A]=get_file_type(filecell{1,1});
     140    FileType_B=FileType_A;
     141    MovieObject_B=MovieObject_A;
     142    if size(filecell,1)>=2 && ~strcmp(filecell{1,1},filecell{2,1})
     143        [FileType_B,FileInfo,MovieObject_B]=get_file_type(filecell{2,1});
     144        CheckImage_B=~isempty(find(strcmp(FileType,ImageTypeOptions)));% =1 for images
     145    end
     146end
     147
     148%%%%% MAIN LOOP %%%%%%
     149
     150MovieObject_A=[];
     151for ifield=1:NbField
     152   
     153    %% Civ1
     154    if isfield (Param.ActionInput,'Civ1')
     155        par_civ1=Param.ActionInput.Civ1;
     156        if isfield(par_civ1,'reverse_pair')% A REVOIR
     157            if par_civ1.reverse_pair
     158                if ischar(par_civ1.ImageB)
     159                    temp=par_civ1.ImageA;
     160                    par_civ1.ImageA=imread(par_civ1.ImageB);
     161                end
     162                if ischar(temp)
     163                    par_civ1.ImageB=imread(temp);
     164                end
     165            end
     166        else
     167            %         if ~isfield(Param.Civ1,'ImageA')
     168            ImageName_A=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,i1_series_Civ1(ifield),[],j1_series_Civ1(ifield));
     169            [par_civ1.ImageA,MovieObject_A] = read_image(ImageName_A,FileType_A,MovieObject_A,FrameIndex_A_Civ1);
     170            %         elseif ischar(Param.Civ1.ImageA)
     171            %             Param.Civ1.ImageA=regexprep(Param.Civ1.ImageA,'''','\');
     172            %             [par_civ1.ImageA,VideoObject] = read_image(Param.Civ1.ImageA,par_civ1.FileTypeA,MovieObject_A,par_civ1.FrameIndexA);
     173            %         end
     174            %         if ~isfield(Param.Civ1,'ImageB')
     175            ImageName_B=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,i2_series_Civ1(ifield),[],j2_series_Civ1(ifield));
     176            [par_civ1.ImageB,MovieObject_B] = read_image(ImageName_B,FileType_B,MovieObject_B,FrameIndex_B_Civ1);
     177            %         elseif isfield(Param.Civ1,'ImageB')&& ischar(Param.Civ1.ImageB)
     178            %              Param.Civ1.ImageB=regexprep(Param.Civ1.ImageB,'''','\');
     179            %              if strcmp(Param.Civ1.ImageA,Param.Civ1.ImageB)% use the same movie object
     180            %                  [par_civ1.ImageB,VideoObject] = read_image(Param.Civ1.ImageB,par_civ1.FileTypeB,VideoObject,par_civ1.FrameIndexB);
     181            %              else
     182            %             [par_civ1.ImageB,VideoObject] = read_image(Param.Civ1.ImageB,par_civ1.FileTypeB,par_civ1.ImageB,par_civ1.FrameIndexB);
     183            %              end
     184            %         end
     185        end
     186        ncfile=fullfile_uvmat(RootPath,OutputDir,RootFile,'.nc',NomTypeNc,i1_series_Civ1(ifield),i2_series_Civ1(ifield),...
     187            j1_series_Civ1(ifield),j2_series_Civ1(ifield));
     188        par_civ1.ImageWidth=FileInfo.Width;
     189        par_civ1.ImageHeight=FileInfo.Height;
     190        list_param=(fieldnames(Param.ActionInput.Civ1))';
     191        Civ1_param=list_param;%default
    193192       
    194         %% transform the input field (e.g; phys) if requested
    195         if ~isempty(transform_fct)
    196             if nargin(transform_fct)>=2
    197                 Data{iview}=transform_fct(Data{iview},XmlData{iview});
    198             else
    199                 Data{iview}=transform_fct(Data{iview});
     193        %set the values of all the global attributes in list_param
     194        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param];
     195        for ilist=1:length(list_param)
     196            Civ1_param{ilist}=['Civ1_' list_param{ilist}];
     197            Data.(['Civ1_' list_param{ilist}])=Param.ActionInput.Civ1.(list_param{ilist});
     198        end
     199        Data.CivStage=1;
     200       
     201        % set the list of variables
     202        Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_F','Civ1_C'};%  cell array containing the names of the fields to record
     203        Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
     204        Data.VarAttribute{1}.Role='coord_x';
     205        Data.VarAttribute{2}.Role='coord_y';
     206        Data.VarAttribute{3}.Role='vector_x';
     207        Data.VarAttribute{4}.Role='vector_y';
     208        Data.VarAttribute{5}.Role='warnflag';
     209       
     210        if strcmp(Param.ActionInput.ListCompareMode, 'PIV volume')
     211            Data.ListVarName=[Data.ListVarName 'Civ1_Z'];
     212            Data.Civ1_X=[];Data.Civ1_Y=[];Data.Civ1_Z=[];
     213            Data.Civ1_U=[];Data.Civ1_V=[];Data.Civ1_C=[];Data.Civ1_F=[];
     214            for ivol=1:NbSlice
     215                % caluclate velocity data (y and v in indices, reverse to y component)
     216                [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
     217                if ~isempty(errormsg)
     218                    return
     219                end
     220                Data.Civ1_X=[Data.Civ1_X reshape(xtable,[],1)];
     221                Data.Civ1_Y=[Data.Civ1_Y reshape(Param.Civ1.ImageHeight-ytable+1,[],1)];
     222                Data.Civ1_Z=[Data.Civ1_Z ivol*ones(numel(xtable),1)];% z=image index in image coordinates
     223                Data.Civ1_U=[Data.Civ1_U reshape(utable,[],1)];
     224                Data.Civ1_V=[Data.Civ1_V reshape(-vtable,[],1)];
     225                Data.Civ1_C=[Data.Civ1_C reshape(ctable,[],1)];
     226                Data.Civ1_F=[Data.Civ1_C reshape(F,[],1)];
    200227            end
    201         end
    202        
    203         %% check whether tps is needed, then calculate tps coefficients if needed
    204         check_tps=0;
    205         if isfield(Param.InputFields,'FieldName')
    206             if ischar(Param.InputFields.FieldName)
    207                 Param.InputFields.FieldName={Param.InputFields.FieldName};
    208             end
    209         else
    210             Param.InputFields.FieldName={};
    211         end
    212         for ilist=1:numel(Param.InputFields.FieldName)
    213             switch Param.InputFields.FieldName{ilist}
    214                 case {'vort','div','strain'}
    215                     check_tps=1;
    216             end
    217         end
    218        
    219         %% calculate tps coeff if needed
    220         check_proj_tps= ~isempty(Param.ProjObject)&& strcmp(Param.ProjObject.ProjMode,'filter')&&~isfield(Data{iview},'Coord_tps');
    221         Data{iview}=tps_coeff_field(Data{iview},check_proj_tps);
    222        
    223         %% projection on object (gridded plane)
    224         if Param.CheckObject
    225             [Data{iview},errormsg]=proj_field(Data{iview},Param.ProjObject);
     228        else %usual PIV
     229            % caluclate velocity data (y and v in indices, reverse to y component)
     230            [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
    226231            if ~isempty(errormsg)
    227                 displ_uvmat('ERROR',['error in merge_proge/proj_field: ' errormsg],checkrun)
    228232                return
    229233            end
    230         end
    231     end
    232     %----------END LOOP ON VIEWS----------------------
     234            Data.Civ1_X=reshape(xtable,[],1);
     235            Data.Civ1_Y=reshape(par_civ1.ImageHeight-ytable+1,[],1);
     236            Data.Civ1_U=reshape(utable,[],1);
     237            Data.Civ1_V=reshape(-vtable,[],1);
     238            Data.Civ1_C=reshape(ctable,[],1);
     239            Data.Civ1_F=reshape(F,[],1);
     240        end
     241    else
     242        if exist('ncfile','var')
     243            CivFile=ncfile;
     244        elseif isfield(Param.Patch1,'CivFile')
     245            CivFile=Param.Patch1.CivFile;
     246        end
     247        Data=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
     248        if isfield(Data,'Txt')
     249            errormsg=Data.Txt;
     250            return
     251        end
     252        if ~isempty(Data.absolut_time_T0')%read civx file
     253            check_civx=1;% test for old civx data format
     254            [Data,vardetect,ichoice]=nc2struct(CivFile);%read the variables in the netcdf file
     255        else
     256            Data=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
     257        end
     258    end
    233259   
    234     %% merge the nbview fields
    235     MergeData=merge_field(Data);
    236     if isfield(MergeData,'Txt')
    237         displ_uvmat('ERROR',MergeData.Txt,checkrun)
     260    %% Fix1
     261    if isfield (Param.ActionInput,'Fix1')
     262        ListFixParam=fieldnames(Param.ActionInput.Fix1);
     263        for ilist=1:length(ListFixParam)
     264            ParamName=ListFixParam{ilist};
     265            ListName=['Fix1_' ParamName];
     266            eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
     267            eval(['Data.' ListName '=Param.ActionInput.Fix1.' ParamName ';'])
     268        end
     269        if check_civx
     270            if ~isfield(Data,'fix')
     271                Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix'];
     272                Data.fix=1;
     273                Data.ListVarName=[Data.ListVarName {'vec_FixFlag'}];
     274                Data.VarDimName=[Data.VarDimName {'nb_vectors'}];
     275            end
     276            Data.vec_FixFlag=fix(Param.ActionInput.Fix1,Data.vec_F,Data.vec_C,Data.vec_U,Data.vec_V,Data.vec_X,Data.vec_Y);
     277        else
     278            Data.ListVarName=[Data.ListVarName {'Civ1_FF'}];
     279            Data.VarDimName=[Data.VarDimName {'nb_vec_1'}];
     280            nbvar=length(Data.ListVarName);
     281            Data.VarAttribute{nbvar}.Role='errorflag';
     282            Data.Civ1_FF=fix(Param.ActionInput.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
     283            Data.CivStage=2;
     284        end
     285    end
     286    %% Patch1
     287    if isfield (Param.ActionInput,'Patch1')
     288        if check_civx
     289            errormsg='Civ Matlab input needed for patch';
     290            return
     291        end
     292       
     293        Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
     294        Data.Patch1_FieldSmooth=Param.ActionInput.Patch1.FieldSmooth;
     295        Data.Patch1_MaxDiff=Param.ActionInput.Patch1.MaxDiff;
     296        Data.Patch1_SubDomainSize=Param.ActionInput.Patch1.SubDomainSize;
     297        nbvar=length(Data.ListVarName);
     298        Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
     299        Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
     300            {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
     301        Data.VarAttribute{nbvar+1}.Role='vector_x';
     302        Data.VarAttribute{nbvar+2}.Role='vector_y';
     303        Data.VarAttribute{nbvar+5}.Role='coord_tps';
     304        Data.VarAttribute{nbvar+6}.Role='vector_x';
     305        Data.VarAttribute{nbvar+7}.Role='vector_y';
     306        Data.Civ1_U_smooth=zeros(size(Data.Civ1_X));
     307        Data.Civ1_V_smooth=zeros(size(Data.Civ1_X));
     308        if isfield(Data,'Civ1_FF')
     309            ind_good=find(Data.Civ1_FF==0);
     310        else
     311            ind_good=1:numel(Data.Civ1_X);
     312        end
     313        [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
     314            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);
     315        Data.Civ1_U_smooth(ind_good)=Ures;
     316        Data.Civ1_V_smooth(ind_good)=Vres;
     317        Data.Civ1_FF(ind_good)=FFres;
     318        Data.CivStage=3;
     319    end
     320   
     321    %% Civ2
     322    if isfield (Param.ActionInput,'Civ2')
     323        par_civ2=Param.ActionInput.Civ2;
     324        par_civ2.ImageA=[];
     325        par_civ2.ImageB=[];
     326        %         if ~isfield(Param.Civ1,'ImageA')
     327        ImageName_A_Civ2=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,i1_series_Civ2(ifield),[],j1_series_Civ2(ifield));
     328
     329        if strcmp(ImageName_A_Civ2,ImageName_A) && isequal(FrameIndex_A_Civ1,FrameIndex_A_Civ2)
     330            par_civ2.ImageA=par_civ1.ImageA;
     331        else
     332            [par_civ2.ImageA,MovieObject_A] = read_image(ImageName_A,FileType_A,MovieObject_A,FrameIndex_A_Civ2);
     333        end
     334        ImageName_B_Civ2=fullfile_uvmat(RootPath,SubDir,RootFile,FileExt,NomType,i2_series_Civ2(ifield),[],j2_series_Civ2(ifield));
     335        if strcmp(ImageName_B_Civ2,ImageName_B) && isequal(FrameIndex_B_Civ1,FrameIndex_B_Civ2)
     336            par_civ2.ImageB=par_civ1.ImageB;
     337        else
     338            [par_civ2.ImageB,MovieObject_B] = read_image(ImageName_B,FileType_B,MovieObject_B,FrameIndex_B_Civ2);
     339        end     
     340       
     341        ncfile=fullfile_uvmat(RootPath,OutputDir,RootFile,'.nc',NomTypeNc,i1_series_Civ2(ifield),i2_series_Civ2(ifield),...
     342            j1_series_Civ2(ifield),j2_series_Civ2(ifield));
     343        par_civ2.ImageWidth=FileInfo.Width;
     344        par_civ2.ImageHeight=FileInfo.Height;
     345       
     346        if isfield(par_civ2,'Grid')% grid points set as input file
     347            if ischar(par_civ2.Grid)%read the grid file if the input is a file name
     348                par_civ2.Grid=dlmread(par_civ2.Grid);
     349                par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     350            end
     351        else% automatic grid
     352            minix=floor(par_civ2.Dx/2)-0.5;
     353            maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
     354            miniy=floor(par_civ2.Dy/2)-0.5;
     355            maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
     356            [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
     357            par_civ2.Grid(:,1)=reshape(GridX,[],1);
     358            par_civ2.Grid(:,2)=reshape(GridY,[],1);
     359        end
     360        Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
     361        Shifty=zeros(size(par_civ2.Grid,1),1);
     362        nbval=zeros(size(par_civ2.Grid,1),1);
     363        if par_civ2.CheckDeformation
     364            DUDX=zeros(size(par_civ2.Grid,1),1);
     365            DUDY=zeros(size(par_civ2.Grid,1),1);
     366            DVDX=zeros(size(par_civ2.Grid,1),1);
     367            DVDY=zeros(size(par_civ2.Grid,1),1);
     368        end
     369        NbSubDomain=size(Data.Civ1_SubRange,3);
     370        % get the guess from patch1
     371        for isub=1:NbSubDomain
     372            nbvec_sub=Data.Civ1_NbCentres(isub);
     373            ind_sel=find(GridX>=Data.Civ1_SubRange(1,1,isub) & GridX<=Data.Civ1_SubRange(1,2,isub) & GridY>=Data.Civ1_SubRange(2,1,isub) & GridY<=Data.Civ1_SubRange(2,2,isub));
     374            epoints = [GridX(ind_sel) GridY(ind_sel)];% coordinates of interpolation sites
     375            ctrs=Data.Civ1_Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
     376            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     377            EM = tps_eval(epoints,ctrs);
     378            Shiftx(ind_sel)=Shiftx(ind_sel)+EM*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
     379            Shifty(ind_sel)=Shifty(ind_sel)+EM*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
     380            if par_civ2.CheckDeformation
     381                [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
     382                DUDX(ind_sel)=DUDX(ind_sel)+EMDX*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
     383                DUDY(ind_sel)=DUDY(ind_sel)+EMDY*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
     384                DVDX(ind_sel)=DVDX(ind_sel)+EMDX*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
     385                DVDY(ind_sel)=DVDY(ind_sel)+EMDY*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
     386            end
     387        end
     388        mask='';
     389        if par_civ2.CheckMask&&~isempty(par_civ2.Mask)&& ~strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1
     390            mask=imread(par_civ2.Mask);
     391        end
     392        ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
     393        iby2=ceil(par_civ2.CorrBoxSize(2)/2);
     394        %     isx2=ibx2+4;% search ara +-4 pixels around the guess
     395        %     isy2=iby2+4;
     396        par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
     397        par_civ2.SearchBoxSize(2)=2*iby2+9;
     398        %par_civ2.SearchBoxSize(1)=2*isx2+1;
     399        %par_civ2.SearchBoxSize(2)=2*isy2+1;
     400        par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
     401        par_civ2.Grid=[GridX(nbval>=1)-par_civ2.SearchBoxShift(:,1)/2 GridY(nbval>=1)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
     402        if par_civ2.CheckDeformation
     403            par_civ2.DUDX=DUDX./nbval;
     404            par_civ2.DUDY=DUDY./nbval;
     405            par_civ2.DVDX=DVDX./nbval;
     406            par_civ2.DVDY=DVDY./nbval;
     407        end
     408        % caluclate velocity data (y and v in indices, reverse to y component)
     409        [xtable ytable utable vtable ctable F] = civ (par_civ2);
     410
     411        list_param=(fieldnames(Param.ActionInput.Civ2))';
     412        list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
     413        for ilist=1:length(list_remove)
     414            index=strcmp(list_remove{ilist},list_param);
     415            if ~isempty(find(index,1))
     416                list_param(index)=[];
     417            end
     418        end
     419        for ilist=1:length(list_param)
     420            Civ2_param{ilist}=['Civ2_' list_param{ilist}];
     421            eval(['Data.Civ2_' list_param{ilist} '=Param.ActionInput.Civ2.' list_param{ilist} ';'])
     422        end
     423        if isfield(Data,'Civ2_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile')
     424            Data.Civ1_gridname='';
     425        end
     426        if isfield(Data,'Civ2_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile')
     427            Data.Civ2_maskname='';
     428        end
     429        Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param {'Civ2_Time','Civ2_Dt'}];
     430        nbvar=numel(Data.ListVarName);
     431        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
     432        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
     433        Data.VarAttribute{nbvar+1}.Role='coord_x';
     434        Data.VarAttribute{nbvar+2}.Role='coord_y';
     435        Data.VarAttribute{nbvar+3}.Role='vector_x';
     436        Data.VarAttribute{nbvar+4}.Role='vector_y';
     437        Data.VarAttribute{nbvar+5}.Role='warnflag';
     438        Data.Civ2_X=reshape(xtable,[],1);
     439        Data.Civ2_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
     440        Data.Civ2_U=reshape(utable,[],1);
     441        Data.Civ2_V=reshape(-vtable,[],1);
     442        Data.Civ2_C=reshape(ctable,[],1);
     443        Data.Civ2_F=reshape(F,[],1);
     444        Data.CivStage=Data.CivStage+1;
     445    end
     446   
     447    %% Fix2
     448    if isfield (Param.ActionInput,'Fix2')
     449        ListFixParam=fieldnames(Param.ActionInput.Fix2);
     450        for ilist=1:length(ListFixParam)
     451            ParamName=ListFixParam{ilist};
     452            ListName=['Fix2_' ParamName];
     453            eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
     454            eval(['Data.' ListName '=Param.ActionInput.Fix2.' ParamName ';'])
     455        end
     456        if check_civx
     457            if ~isfield(Data,'fix2')
     458                Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2'];
     459                Data.fix2=1;
     460                Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}];
     461                Data.VarDimName=[Data.VarDimName {'nb_vectors2'}];
     462            end
     463            Data.vec_FixFlag=fix(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y);
     464        else
     465            Data.ListVarName=[Data.ListVarName {'Civ2_FF'}];
     466            Data.VarDimName=[Data.VarDimName {'nb_vec_2'}];
     467            nbvar=length(Data.ListVarName);
     468            Data.VarAttribute{nbvar}.Role='errorflag';
     469            Data.Civ2_FF=fix(Param.ActionInput.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V);
     470            Data.CivStage=Data.CivStage+1;
     471        end
     472       
     473    end
     474   
     475    %% Patch2
     476    if isfield (Param.ActionInput,'Patch2')
     477        Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch2_Rho','Patch2_Threshold','Patch2_SubDomain'}];
     478        Data.Patch2_FieldSmooth=Param.ActionInput.Patch2.FieldSmooth;
     479        Data.Patch2_MaxDiff=Param.ActionInput.Patch2.MaxDiff;
     480        Data.Patch2_SubDomainSize=Param.ActionInput.Patch2.SubDomainSize;
     481        nbvar=length(Data.ListVarName);
     482        Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
     483        Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
     484            {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
     485       
     486        Data.VarAttribute{nbvar+1}.Role='vector_x';
     487        Data.VarAttribute{nbvar+2}.Role='vector_y';
     488        Data.VarAttribute{nbvar+5}.Role='coord_tps';
     489        Data.VarAttribute{nbvar+6}.Role='vector_x';
     490        Data.VarAttribute{nbvar+7}.Role='vector_y';
     491        Data.Civ2_U_smooth=zeros(size(Data.Civ2_X));
     492        Data.Civ2_V_smooth=zeros(size(Data.Civ2_X));
     493        if isfield(Data,'Civ2_FF')
     494            ind_good=find(Data.Civ2_FF==0);
     495        else
     496            ind_good=1:numel(Data.Civ2_X);
     497        end
     498        [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures, Vres,tild,FFres]=...
     499            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);
     500        Data.Civ2_U_smooth(ind_good)=Ures;
     501        Data.Civ2_V_smooth(ind_good)=Vres;
     502        Data.Civ2_FF(ind_good)=FFres;
     503        Data.CivStage=Data.CivStage+1;
     504    end
     505   
     506    %% write result in a netcdf file if requested
     507    if exist('ncfile','var')
     508        errormsg=struct2nc(ncfile,Data);
     509        if isempty(errormsg)
     510            disp([ncfile ' written'])
     511        else
     512            disp(errormsg)
     513        end
     514    end
     515end
     516
     517% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
     518%--------------------------------------------------------------------------
     519% function [xtable ytable utable vtable typevector] = civ (image1,image2,ibx,iby step, subpixfinder, mask, roi)
     520%
     521% OUTPUT:
     522% xtable: set of x coordinates
     523% ytable: set of y coordiantes
     524% utable: set of u displacements (along x)
     525% vtable: set of v displacements (along y)
     526% ctable: max image correlation for each vector
     527% typevector: set of flags, =1 for good, =0 for NaN vectors
     528%
     529%INPUT:
     530% par_civ: structure of input parameters, with fields:
     531%  .CorrBoxSize
     532%  .SearchBoxSize
     533%  .SearchBoxShift
     534%  .ImageHeight
     535%  .ImageWidth
     536%  .Dx, Dy
     537%  .Grid
     538%  .Mask
     539%  .MinIma
     540%  .MaxIma
     541%  .image1:first image (matrix)
     542% image2: second image (matrix)
     543% ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1)
     544% isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1)
     545% shiftx, shifty: shift of the search box (in pixel index, yshift reversed)
     546% step: mesh of the measurement points (in px)
     547% subpixfinder=1 or 2 controls the curve fitting of the image correlation
     548% mask: =[] for no mask
     549% roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[];
     550function [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ)
     551%this funtion performs the DCC PIV analysis. Recent window-deformation
     552%methods perform better and will maybe be implemented in the future.
     553
     554%% prepare measurement grid
     555if isfield(par_civ,'Grid')% grid points set as input
     556    if ischar(par_civ.Grid)%read the drid file if the input is a file name
     557        par_civ.Grid=dlmread(par_civ.Grid);
     558        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
     559    end
     560else% automatic grid
     561    minix=floor(par_civ.Dx/2)-0.5;
     562    maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
     563    miniy=floor(par_civ.Dy/2)-0.5;
     564    maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
     565    [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
     566    par_civ.Grid(:,1)=reshape(GridX,[],1);
     567    par_civ.Grid(:,2)=reshape(GridY,[],1);
     568end
     569nbvec=size(par_civ.Grid,1);
     570
     571%% prepare correlation and search boxes
     572ibx2=ceil(par_civ.CorrBoxSize(1)/2);
     573iby2=ceil(par_civ.CorrBoxSize(2)/2);
     574isx2=ceil(par_civ.SearchBoxSize(1)/2);
     575isy2=ceil(par_civ.SearchBoxSize(2)/2);
     576shiftx=round(par_civ.SearchBoxShift(:,1));
     577shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
     578if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
     579    shiftx=shiftx*ones(nbvec,1);
     580    shifty=shifty*ones(nbvec,1);
     581end
     582
     583%% Default output
     584xtable=par_civ.Grid(:,1);
     585ytable=par_civ.Grid(:,2);
     586utable=zeros(nbvec,1);
     587vtable=zeros(nbvec,1);
     588ctable=zeros(nbvec,1);
     589F=zeros(nbvec,1);
     590result_conv=[];
     591errormsg='';
     592
     593%% prepare mask
     594if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
     595    if strcmp(par_civ.Mask,'all')
     596        return    % get the grid only, no civ calculation
     597    elseif ischar(par_civ.Mask)
     598        par_civ.Mask=imread(par_civ.Mask);
     599    end
     600end
     601check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
     602check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
     603
     604% %% prepare images
     605% if isfield(par_civ,'reverse_pair')
     606%     if par_civ.reverse_pair
     607%         if ischar(par_civ.ImageB)
     608%             temp=par_civ.ImageA;
     609%             par_civ.ImageA=imread(par_civ.ImageB);
     610%         end
     611%         if ischar(temp)
     612%             par_civ.ImageB=imread(temp);
     613%         end
     614%     end
     615% else
     616%     if ischar(par_civ.ImageA)
     617%         par_civ.ImageA=imread(par_civ.ImageA);
     618%     end
     619%     if ischar(par_civ.ImageB)
     620%         par_civ.ImageB=imread(par_civ.ImageB);
     621%     end
     622% end
     623par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
     624par_civ.ImageB=sum(double(par_civ.ImageB),3);
     625[npy_ima npx_ima]=size(par_civ.ImageA);
     626if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima])
     627    errormsg='image pair with unequal size';
     628    return
     629end
     630
     631%% Apply mask
     632    % Convention for mask IDEAS TO IMPLEMENT ?
     633    % mask >200 : velocity calculated
     634    %  200 >=mask>150;velocity not calculated, interpolation allowed (bad spots)
     635    % 150>=mask >100: velocity not calculated, nor interpolated
     636    %  100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries)
     637    %  20>=mask: velocity=0
     638checkmask=0;
     639MinA=min(min(par_civ.ImageA));
     640MinB=min(min(par_civ.ImageB));
     641if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
     642   checkmask=1;
     643   if ~isequal(size(par_civ.Mask),[npy_ima npx_ima])
     644        errormsg='mask must be an image with the same size as the images';
    238645        return
    239     end
     646   end
     647  %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
     648    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
     649    par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined  area
     650    par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
     651end
     652
     653%% compute image correlations: MAINLOOP on velocity vectors
     654corrmax=0;
     655sum_square=1;% default
     656mesh=1;% default
     657CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1;
     658if CheckDecimal
     659    mesh=0.2;%mesh in pixels for subpixel image interpolation
     660    CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
     661end
     662% vector=[0 0];%default
     663for ivec=1:nbvec
     664    iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
     665    jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
     666    %if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
     667    %         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
     668    %               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
     669    %             F(ivec)=3;
     670    %         else
     671    F(ivec)=0;
     672    subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
     673    subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
     674    subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
     675    subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
     676    image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
     677    image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
     678    check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
     679    check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
     680    check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
     681    check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
    240682   
    241     % time of the merged field:
    242     if ~isempty(time)% time defined from ImaDoc
    243         timeread=time(:,index);
    244     end
    245     timeread=mean(timeread);
    246    
    247     % generating the name of the merged field
    248     i1=i1_series{iview}(index);
    249     if ~isempty(i2_series{iview})
    250         i2=i2_series{iview}(index);
     683    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
     684    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
     685    image1_mean=mean(mean(image1_crop));
     686    image2_mean=mean(mean(image2_crop));
     687    %threshold on image minimum
     688    if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
     689        F(ivec)=3;
     690    end
     691    %threshold on image maximum
     692    if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
     693        F(ivec)=3;
     694    end
     695    %         end
     696    if F(ivec)~=3
     697        image1_crop=image1_crop-image1_mean;%substract the mean
     698        image2_crop=image2_crop-image2_mean;
     699        if CheckDecimal
     700            xi=(1:mesh:size(image1_crop,2));
     701            yi=(1:mesh:size(image1_crop,1))';
     702            if CheckDeformation
     703                [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
     704                XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
     705                YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
     706                image1_crop=interp2(image1_crop,XIant,YIant);
     707            else
     708                image1_crop=interp2(image1_crop,xi,yi);
     709            end
     710            xi=(1:mesh:size(image2_crop,2));
     711            yi=(1:mesh:size(image2_crop,1))';
     712            image2_crop=interp2(image2_crop,xi,yi);
     713        end
     714        sum_square=sum(sum(image1_crop.*image1_crop));
     715        %reference: Oliver Pust, PIV: Direct Cross-Correlation
     716        result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
     717        corrmax= max(max(result_conv));
     718        result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
     719        %Find the correlation max, at 255
     720        [y,x] = find(result_conv==255,1);
     721        if ~isempty(y) && ~isempty(x)
     722            try
     723                if par_civ.CorrSmooth==1
     724                    [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
     725                elseif par_civ.CorrSmooth==2
     726                    [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
     727                end
     728                utable(ivec)=vector(1)*mesh+shiftx(ivec);
     729                vtable(ivec)=vector(2)*mesh+shifty(ivec);
     730                xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
     731                ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
     732                iref=round(xtable(ivec));% image index for the middle of the vector
     733                jref=round(ytable(ivec));
     734                if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
     735                    utable(ivec)=0;
     736                    vtable(ivec)=0;
     737                    F(ivec)=3;
     738                end
     739                ctable(ivec)=corrmax/sum_square;% correlation value
     740            catch ME
     741                F(ivec)=3;
     742            end
     743        else
     744            F(ivec)=3;
     745        end
     746    end
     747end
     748result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
     749
     750%------------------------------------------------------------------------
     751% --- Find the maximum of the correlation function after interpolation
     752function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
     753%------------------------------------------------------------------------
     754vector=[0 0]; %default
     755F=0;
     756[npy,npx]=size(result_conv);
     757
     758% if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
     759    %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ï¿œ Israel Institute of Technology
     760    %http://urapiv.wordpress.com
     761    peaky = y;
     762    if y <= npy-1 && y >= 1
     763        f0 = log(result_conv(y,x));
     764        f1 = real(log(result_conv(y-1,x)));
     765        f2 = real(log(result_conv(y+1,x)));
     766        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
    251767    else
    252         i2=i1;
    253     end
    254     j1=1;
    255     j2=1;
    256     if ~isempty(j1_series{iview})
    257         j1=j1_series{iview}(index);
    258         if ~isempty(j2_series{iview})
    259             j2=j2_series{iview}(index);
    260         else
    261             j2=j1;
    262         end
    263     end
    264     OutputFile=fullfile_uvmat(RootPath{1},OutputSubDir,RootFile{1},FileExtOut,NomType{1},i1,i2,j1,j2);
    265    
    266     % recording the merged field
    267     if CheckImage{1}    %in case of input images an image is produced
    268         if isa(MergeData.A,'uint8')
    269             bitdepth=8;
    270         elseif isa(MergeData.A,'uint16')
    271             bitdepth=16;
    272         end
    273         imwrite(MergeData.A,OutputFile,'BitDepth',bitdepth);
    274         %write xml calibration file
    275         siz=size(MergeData.A);
    276         npy=siz(1);
    277         npx=siz(2);
    278         if isfield(MergeData,'VarAttribute')&&isfield(MergeData.VarAttribute{1},'Coord_2')&&isfield(MergeData.VarAttribute{1},'Coord_1')
    279             Rangx=MergeData.VarAttribute{1}.Coord_2;
    280             Rangy=MergeData.VarAttribute{1}.Coord_1;
    281         elseif isfield(MergeData,'AX')&& isfield(MergeData,'AY')
    282             Rangx=[MergeData.AX(1) MergeData.AX(end)];
    283             Rangy=[MergeData.AY(1) MergeData.AY(end)];
    284         else
    285             Rangx=[0.5 npx-0.5];
    286             Rangy=[npy-0.5 0.5];%default
    287         end
    288         pxcmx=(npx-1)/(Rangx(2)-Rangx(1));
    289         pxcmy=(npy-1)/(Rangy(1)-Rangy(2));
    290         T_x=-pxcmx*Rangx(1)+0.5;
    291         T_y=-pxcmy*Rangy(2)+0.5;
    292         GeometryCal.focal=1;
    293         GeometryCal.R=[pxcmx,0,0;0,pxcmy,0;0,0,1];
    294         GeometryCal.Tx_Ty_Tz=[T_x T_y 1];
    295         ImaDoc.GeometryCalib=GeometryCal;
    296         %             t=struct2xml(ImaDoc);
    297         %             t=set(t,1,'name','ImaDoc');
    298         %             save(t,[filebase_merge '.xml'])
    299         %             display([filebase_merge '.xml saved'])
     768        F=-2; % warning flag for vector truncated by the limited search box
     769    end
     770    peakx=x;
     771    if x <= npx-1 && x >= 1
     772        f0 = log(result_conv(y,x));
     773        f1 = real(log(result_conv(y,x-1)));
     774        f2 = real(log(result_conv(y,x+1)));
     775        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
    300776    else
    301         MergeData.ListGlobalAttribute={'Conventions','Project','InputFile_1','InputFile_end','nb_coord','nb_dim','dt','Time','civ'};
    302         MergeData.Conventions='uvmat';
    303         MergeData.nb_coord=2;
    304         MergeData.nb_dim=2;
    305         dt=[];
    306         if isfield(Data{1},'dt')&& isnumeric(Data{1}.dt)
    307             dt=Data{1}.dt;
    308         end
    309         for iview =2:numel(Data)
    310             if ~(isfield(Data{iview},'dt')&& isequal(Data{iview}.dt,dt))
    311                 dt=[];%dt not the same for all fields
    312             end
    313         end
    314         if isempty(dt)
    315             MergeData.ListGlobalAttribute(6)=[];
    316         else
    317             MergeData.dt=dt;
    318         end
    319         MergeData.Time=timeread;
    320         error=struct2nc(OutputFile,MergeData);%save result file
    321         if isempty(error)
    322             display(['output file ' OutputFile ' written'])
    323         else
    324             display(error)
    325         end
    326     end
    327 end
    328 
    329 %'merge_field': concatene fields
     777        F=-2; % warning flag for vector truncated by the limited search box
     778    end
     779    vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     780
    330781%------------------------------------------------------------------------
    331 function MergeData=merge_field(Data)
    332 %% default output
    333 if isempty(Data)||~iscell(Data)
    334     MergeData=[];
    335     return
    336 end
    337 MergeData=Data{1};%default
    338 error=0;
    339 nbview=length(Data);
    340 if nbview==1
    341     return
    342 end
    343 
    344 %% group the variables (fields of 'Data') in cells of variables with the same dimensions
    345 [CellVarIndex,NbDim,VarTypeCell]=find_field_cells(Data{1});
    346 %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
    347 % CellVarIndex=cells of variable index arrays
    348 for icell=1:length(CellVarIndex)
    349     if NbDim(icell)==1
    350         continue
    351     end
    352     VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
    353     VarType=VarTypeCell{icell};
    354     ivar_X=VarType.coord_x;
    355     ivar_Y=VarType.coord_y;
    356     ivar_FF=VarType.errorflag;
    357     if isempty(ivar_X)
    358         test_grid=1;%test for input data on regular grid (e.g. image)coordinates
     782% --- Find the maximum of the correlation function after interpolation
     783function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
     784%------------------------------------------------------------------------
     785vector=[0 0]; %default
     786F=-2;
     787peaky=y;
     788peakx=x;
     789[npy,npx]=size(result_conv);
     790if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1)
     791    F=0;
     792    for i=-1:1
     793        for j=-1:1
     794            %following 15 lines based on
     795            %H. Nobach ï¿œ M. Honkanen (2005)
     796            %Two-dimensional Gaussian regression for sub-pixel displacement
     797            %estimation in particle image velocimetry or particle position
     798            %estimation in particle tracking velocimetry
     799            %Experiments in Fluids (2005) 38: 511ï¿œ515
     800            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
     801            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
     802            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
     803            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
     804            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
     805        end
     806    end
     807    c10=(1/6)*sum(sum(c10));
     808    c01=(1/6)*sum(sum(c01));
     809    c11=(1/4)*sum(sum(c11));
     810    c20=(1/6)*sum(sum(c20));
     811    c02=(1/6)*sum(sum(c02));
     812    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
     813    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
     814    if abs(deltax)<1
     815        peakx=x+deltax;
     816    end
     817    if abs(deltay)<1
     818        peaky=y+deltay;
     819    end
     820end
     821vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
     822
     823%'RUN_FIX': function for fixing velocity fields:
     824%-----------------------------------------------
     825% RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)
     826%
     827%filename: name of the netcdf file (used as input and output)
     828%field: structure specifying the names of the fields to fix (depending on civ1 or civ2)
     829    %.vel_type='civ1' or 'civ2';
     830    %.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1);
     831    %.fixflag=name of fix flag variable ('vec_FixFlag' for civ1)
     832%flagindex: flag specifying which values of vec_f are removed:
     833        % if flagindex(1)=1: vec_f=-2 vectors are removed
     834        % if flagindex(2)=1: vec_f=3 vectors are removed
     835        % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
     836%iter=1 for civ1 fields and iter=2 for civ2 fields
     837%thresh_vecC: threshold in the image correlation vec_C
     838%flag_mask: =1 mask used to remove vectors (0 else)
     839%maskname: name of the mask image file for fix
     840%thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
     841%inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
     842%fileref: .nc file name for a reference velocity (='': refrence 0 used)
     843%fieldref: 'civ1','filter1'...feld used in fileref
     844
     845function FF=fix(Param,F,C,U,V,X,Y)
     846FF=zeros(size(F));%default
     847
     848%criterium on warn flags
     849FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'};
     850FlagVal=[-2 2 3 4];
     851for iflag=1:numel(FlagName)
     852    if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag})
     853        FF=(FF==1| F==FlagVal(iflag));
     854    end
     855end
     856%criterium on correlation values
     857if isfield (Param,'MinCorr')
     858    FF=FF==1 | C<Param.MinCorr;
     859end
     860if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
     861    Umod= U.*U+V.*V;
     862    if isfield (Param,'MinVel')&&~isempty(Param.MinVel)
     863        FF=FF==1 | Umod<(Param.MinVel*Param.MinVel);
     864    end
     865    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
     866        FF=FF==1 | Umod>(Param.MaxVel*Param.MaxVel);
     867    end
     868end
     869
     870
     871%------------------------------------------------------------------------
     872% --- determine the list of index pairs of processing file
     873function [i1_series,i2_series,j1_series,j2_series,check_bounds,NomTypeNc]=...
     874    find_pair_indices(str_civ,i_series,j_series,MinIndex,MaxIndex)
     875%------------------------------------------------------------------------
     876i1_series=i_series;% set of first image numbers
     877i2_series=i_series;
     878j1_series=ones(size(i_series));% set of first image numbers
     879j2_series=ones(size(i_series));
     880check_bounds=false(size(i_series));
     881r=regexp(str_civ,'^\D(?<ind>[i|j])= (?<num1>\d+)\|(?<num2>\d+)','names');
     882if ~isempty(r)
     883    mode=['D' r.ind];
     884    ind1=str2num(r.num1);
     885    ind2=str2num(r.num2);
     886else
     887    mode='burst';
     888    r=regexp(str_civ,'^j= (?<num1>[a-z])-(?<num2>[a-z])','names');
     889    if ~isempty(r)
     890        NomTypeNc='_1ab';
    359891    else
    360         if length(ivar_Y)~=1
    361                 displ_uvmat('ERROR','y coordinate missing in proj_field.m',checkrun)
    362                 return
    363         end
    364         test_grid=0;
    365     end
    366     %case of input fields with unstructured coordinates
    367     if ~test_grid
    368         for ivar=VarIndex
    369             VarName=MergeData.ListVarName{ivar};
    370             for iview=1:nbview
    371                 MergeData.(VarName)=[MergeData.(VarName); Data{iview}.(VarName)];
    372             end
    373         end
    374     %case of fields defined on a structured  grid
    375     else 
    376         testFF=0;
    377         for iview=2:nbview
    378             for ivar=VarIndex
    379                 VarName=MergeData.ListVarName{ivar};
    380                 if isfield(MergeData,'VarAttribute')
    381                     if length(MergeData.VarAttribute)>=ivar && isfield(MergeData.VarAttribute{ivar},'Role') && isequal(MergeData.VarAttribute{ivar}.Role,'errorflag')
    382                         testFF=1;
    383                     end
    384                 end
    385                 MergeData.(VarName)=MergeData.(VarName) + Data{iview}.(VarName);
    386             end
    387         end
    388         if testFF
    389             nbaver=nbview-MergeData.FF;
    390             indgood=find(nbaver>0);
    391             for ivar=VarIndex
    392                 VarName=MergeData.ListVarName{ivar};
    393                 MergeData.(VarName)(indgood)=double(MergeData.(VarName)(indgood))./nbaver(indgood);
    394             end
    395         else
    396             for ivar=VarIndex
    397                 VarName=MergeData.ListVarName{ivar};
    398                 MergeData.(VarName)=double(MergeData.(VarName))./nbview;
    399             end   
    400         end
    401     end
    402 end
    403 
    404    
     892        r=regexp(str_civ,'^j= (?<num1>[A-Z])-(?<num2>[A-Z])','names');
     893        if ~isempty(r)
     894            NomTypeNc='_1AB';
     895        else
     896            r=regexp(str_civ,'^j= (?<num1>\d+)-(?<num2>\d+)','names');
     897            if ~isempty(r)
     898                NomTypeNc='_1_1-2';
     899            end           
     900        end
     901    end
     902    if isempty(r)
     903        display('wrong pair mode input option')
     904    else
     905    ind1=stra2num(r.num1);
     906    ind2=stra2num(r.num2);
     907    end
     908end
     909if strcmp (mode,'Di')
     910    i1_series=i_series-ind1;% set of first image numbers
     911    i2_series=i_series+ind2;
     912    check_bounds=i1_series<MinIndex(1,1) | i2_series>MaxIndex(1,1);
     913    if isempty(j_series)
     914        NomTypeNc='_1-2';
     915    else
     916        NomTypeNc='_1-2_1';
     917    end
     918elseif strcmp (mode,'Dj')
     919    j1_series=j_series-ind1;
     920    j2_series=j_series+ind2;
     921    check_bounds=j1_series<MinIndex(1,2) | j2_series>MaxIndex(1,2);
     922    NomTypeNc='_1_1-2';
     923else  %bursts
     924    j1_series=ind1*ones(size(i_series));
     925    j2_series=ind2*ones(size(i_series));
     926end
     927
     928%     if length(indsel)>=1
     929%         firstind=indsel(1);
     930%         lastind=indsel(end);
     931%         set(handles.first_j,'String',num2str(ref_j(firstind)))%update the display of first and last fields
     932%         set(handles.last_j,'String',num2str(ref_j(lastind)))
     933%         ref_j=ref_j(indsel);
     934%         j1_civ1=j1_civ1(indsel);
     935%         j2_civ1=j2_civ1(indsel);
     936%         j1_civ2=j1_civ2(indsel);
     937%         j2_civ2=j2_civ2(indsel);
     938%     end
     939% elseif isequal(mode,'pair j1-j2') %case of bursts (png_old or png_2D)
     940%     displ_num=get(handles.ListPairCiv1,'UserData');
     941%     i1_civ1=ref_i;
     942%     i2_civ1=ref_i;
     943%     j1_civ1=displ_num(1,index_civ1);
     944%     j2_civ1=displ_num(2,index_civ1);
     945%     i1_civ2=ref_i;
     946%     i2_civ2=ref_i;
     947%     j1_civ2=displ_num(1,index_civ2);
     948%     j2_civ2=displ_num(2,index_civ2);
     949% elseif isequal(mode,'displacement')
     950%     i1_civ1=ref_i;
     951%     i2_civ1=ref_i;
     952%     j1_civ1=ref_j;
     953%     j2_civ1=ref_j;
     954%     i1_civ2=ref_i;
     955%     i2_civ2=ref_i;
     956%     j1_civ2=ref_j;
     957%     j2_civ2=ref_j;
     958% end
     959
     960
     961
     962
Note: See TracChangeset for help on using the changeset viewer.