Changeset 598 for trunk/src/series/civ_series.m
 Timestamp:
 Apr 2, 2013, 9:13:42 AM (12 years ago)
 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 subfunctions: 3 % civ: PIV function itself 4 % fix: removes false vectors after detection by various criteria 5 % filter_tps: make interpolationsmoothing 2 6 % 3 % function ParamOut=merge_proj(Param) 4 % 5 %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%% 7 % function [Data,errormsg,result_conv]= civ_series(Param,ncfile) 6 8 % 7 9 %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 intercorrelation function for the last grid point (used for tests) 9 13 % 10 14 %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) 15 23 % 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 / CNRSUJFINPG, joel.sommeria@legi.grenobleinp.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 40 function [Data,errormsg,result_conv]= civ_series(Param,ncfile) 41 errormsg=''; 42 path_series=fileparts(which('series')); 43 addpath(fullfile(path_series,'series')) 43 44 %% set the input elements needed on the GUI series when the action is selected in the menu ActionName 44 45 if 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 55 58 end 56 59 … … 63 66 checkrun=0; 64 67 end 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}; 68 OutputDir=[Param.OutputSubDir Param.OutputDirExt]; 69 70 Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; 71 Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes 72 Data.Program='civ_series'; 73 Data.CivStage=0;%default 74 ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'}; %variables to read 75 ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'}; 76 mask=''; 77 maskname='';%default 78 check_civx=0;%default 79 check_civ1=0;%default 80 check_patch1=0;%default 81 82 % case of input Param set by an xml file (batch mode) 83 if ischar(Param) 84 Param=xml2struct(Param); %if Param is the name of an xml file, read this file as a Matlab structure 85 end 86 87 RootPath=Param.InputTable{1,1}; 88 RootFile=Param.InputTable{1,3}; 89 SubDir=Param.InputTable{1,2}; 90 NomType=Param.InputTable{1,4}; 91 FileExt=Param.InputTable{1,5}; 92 PairCiv1=Param.ActionInput.PairIndices.ListPairCiv1; 93 PairCiv2=''; 94 if isfield(Param.ActionInput.PairIndices,'ListPairCiv2') 95 PairCiv2=Param.ActionInput.PairIndices.ListPairCiv2; 96 end 97 98 % option use with GUI series 99 NbField=1; 100 MovieObject_A=[]; 101 if 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; 111 119 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 146 end 147 148 %%%%% MAIN LOOP %%%%%% 149 150 MovieObject_A=[]; 151 for 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 193 192 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.ImageHeightytable+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)]; 200 227 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); 226 231 if ~isempty(errormsg) 227 displ_uvmat('ERROR',['error in merge_proge/proj_field: ' errormsg],checkrun)228 232 return 229 233 end 230 end 231 end 232 %END LOOP ON VIEWS 234 Data.Civ1_X=reshape(xtable,[],1); 235 Data.Civ1_Y=reshape(par_civ1.ImageHeightytable+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 233 259 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.ImageWidth1)/par_civ2.Dx); 354 miniy=floor(par_civ2.Dy/2)0.5; 355 maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight1)/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 515 end 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=[]; 550 function [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ) 551 %this funtion performs the DCC PIV analysis. Recent windowdeformation 552 %methods perform better and will maybe be implemented in the future. 553 554 %% prepare measurement grid 555 if 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 560 else% automatic grid 561 minix=floor(par_civ.Dx/2)0.5; 562 maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth1)/par_civ.Dx); 563 miniy=floor(par_civ.Dy/2)0.5; 564 maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight1)/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); 568 end 569 nbvec=size(par_civ.Grid,1); 570 571 %% prepare correlation and search boxes 572 ibx2=ceil(par_civ.CorrBoxSize(1)/2); 573 iby2=ceil(par_civ.CorrBoxSize(2)/2); 574 isx2=ceil(par_civ.SearchBoxSize(1)/2); 575 isy2=ceil(par_civ.SearchBoxSize(2)/2); 576 shiftx=round(par_civ.SearchBoxShift(:,1)); 577 shifty=round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases 578 if 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); 581 end 582 583 %% Default output 584 xtable=par_civ.Grid(:,1); 585 ytable=par_civ.Grid(:,2); 586 utable=zeros(nbvec,1); 587 vtable=zeros(nbvec,1); 588 ctable=zeros(nbvec,1); 589 F=zeros(nbvec,1); 590 result_conv=[]; 591 errormsg=''; 592 593 %% prepare mask 594 if 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 600 end 601 check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold 602 check_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 623 par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images 624 par_civ.ImageB=sum(double(par_civ.ImageB),3); 625 [npy_ima npx_ima]=size(par_civ.ImageA); 626 if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima]) 627 errormsg='image pair with unequal size'; 628 return 629 end 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 638 checkmask=0; 639 MinA=min(min(par_civ.ImageA)); 640 MinB=min(min(par_civ.ImageB)); 641 if 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'; 238 645 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 651 end 652 653 %% compute image correlations: MAINLOOP on velocity vectors 654 corrmax=0; 655 sum_square=1;% default 656 mesh=1;% default 657 CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1; 658 if CheckDecimal 659 mesh=0.2;%mesh in pixels for subpixel image interpolation 660 CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1; 661 end 662 % vector=[0 0];%default 663 for 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.ImageHeightpar_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 jrefiby2<1  jref+iby2>par_civ.ImageHeight irefibx2<1  iref+ibx2>par_civ.ImageWidth... 668 % jref+shifty(ivec)isy2<1jref+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=irefibx2:iref+ibx2;% x indices defining the first subimage 673 subrange1_y=jrefiby2: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; 240 682 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_cropimage1_mean;%substract the mean 698 image2_crop=image2_cropimage2_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(xiceil(size(image1_crop,2)/2),yiceil(size(image1_crop,1)/2)); 704 XIant=XIpar_civ.DUDX(ivec)*XIpar_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2); 705 YIant=YIpar_civ.DVDX(ivec)*XIpar_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 CrossCorrelation 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)/20.5;% convec flow (velocity taken at the point middle from imgae 1 and 2) 731 ytable(ivec)=jref+vtable(ivec)/20.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 747 end 748 result_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 752 function [vector,F] = SUBPIXGAUSS (result_conv,x,y) 753 % 754 vector=[0 0]; %default 755 F=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 <= npy1 && y >= 1 763 f0 = log(result_conv(y,x)); 764 f1 = real(log(result_conv(y1,x))); 765 f2 = real(log(result_conv(y+1,x))); 766 peaky = peaky+ (f1f2)/(2*f14*f0+2*f2); 251 767 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 npx0.5]; 286 Rangy=[npy0.5 0.5];%default 287 end 288 pxcmx=(npx1)/(Rangx(2)Rangx(1)); 289 pxcmy=(npy1)/(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 <= npx1 && x >= 1 772 f0 = log(result_conv(y,x)); 773 f1 = real(log(result_conv(y,x1))); 774 f2 = real(log(result_conv(y,x+1))); 775 peakx = peakx+ (f1f2)/(2*f14*f0+2*f2); 300 776 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=[peakxfloor(npx/2)1 peakyfloor(npy/2)1]; 780 330 781 % 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 783 function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y) 784 % 785 vector=[0 0]; %default 786 F=2; 787 peaky=y; 788 peakx=x; 789 [npy,npx]=size(result_conv); 790 if (x <= npx1) && (y <= npy1) && (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 %Twodimensional Gaussian regression for subpixel 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^22)*log(result_conv(y+j, x+i)); 804 c02(j+2,i+2)=(3*j^22)*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*c012*c10*c02)/(4*c20*c02c11^2); 813 deltay=(c11*c102*c01*c20)/(4*c20*c02c11^2); 814 if abs(deltax)<1 815 peakx=x+deltax; 816 end 817 if abs(deltay)<1 818 peaky=y+deltay; 819 end 820 end 821 vector=[peakxfloor(npx/2)1 peakyfloor(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 845 function FF=fix(Param,F,C,U,V,X,Y) 846 FF=zeros(size(F));%default 847 848 %criterium on warn flags 849 FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'}; 850 FlagVal=[2 2 3 4]; 851 for iflag=1:numel(FlagName) 852 if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag}) 853 FF=(FF==1 F==FlagVal(iflag)); 854 end 855 end 856 %criterium on correlation values 857 if isfield (Param,'MinCorr') 858 FF=FF==1  C<Param.MinCorr; 859 end 860 if (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 868 end 869 870 871 % 872 %  determine the list of index pairs of processing file 873 function [i1_series,i2_series,j1_series,j2_series,check_bounds,NomTypeNc]=... 874 find_pair_indices(str_civ,i_series,j_series,MinIndex,MaxIndex) 875 % 876 i1_series=i_series;% set of first image numbers 877 i2_series=i_series; 878 j1_series=ones(size(i_series));% set of first image numbers 879 j2_series=ones(size(i_series)); 880 check_bounds=false(size(i_series)); 881 r=regexp(str_civ,'^\D(?<ind>[ij])= (?<num1>\d+)\(?<num2>\d+)','names'); 882 if ~isempty(r) 883 mode=['D' r.ind]; 884 ind1=str2num(r.num1); 885 ind2=str2num(r.num2); 886 else 887 mode='burst'; 888 r=regexp(str_civ,'^j= (?<num1>[az])(?<num2>[az])','names'); 889 if ~isempty(r) 890 NomTypeNc='_1ab'; 359 891 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=nbviewMergeData.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>[AZ])(?<num2>[AZ])','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_12'; 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 908 end 909 if strcmp (mode,'Di') 910 i1_series=i_seriesind1;% 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='_12'; 915 else 916 NomTypeNc='_12_1'; 917 end 918 elseif strcmp (mode,'Dj') 919 j1_series=j_seriesind1; 920 j2_series=j_series+ind2; 921 check_bounds=j1_series<MinIndex(1,2)  j2_series>MaxIndex(1,2); 922 NomTypeNc='_1_12'; 923 else %bursts 924 j1_series=ind1*ones(size(i_series)); 925 j2_series=ind2*ones(size(i_series)); 926 end 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 j1j2') %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.