[309] | 1 | %'civ_matlab': Matlab version of the PIV programs CivX |
---|
| 2 | % --- call the sub-functions: |
---|
| 3 | % civ: PIV function itself |
---|
| 4 | % fix: removes false vectors after detection by various criteria |
---|
| 5 | % patch: make interpolation-smoothing |
---|
| 6 | %------------------------------------------------------------------------ |
---|
| 7 | % function [Data,errormsg,result_conv]= civ_uvmat(Param,ncfile) |
---|
| 8 | % |
---|
| 9 | %OUTPUT |
---|
| 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) |
---|
| 13 | % |
---|
| 14 | %INPUT: |
---|
| 15 | % Param: input images and processing parameters |
---|
| 16 | % ncfile: name of a netcdf file to be created for the result (extension .nc) |
---|
| 17 | % |
---|
| 18 | %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA |
---|
| 19 | % Copyright 2011, LEGI / CNRS-UJF-INPG, joel.sommeria@legi.grenoble-inp.fr. |
---|
| 20 | %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA |
---|
| 21 | % This is part of the toolbox UVMAT. |
---|
| 22 | % |
---|
| 23 | % UVMAT is free software; you can redistribute it and/or modify |
---|
| 24 | % it under the terms of the GNU General Public License as published by |
---|
| 25 | % the Free Software Foundation; either version 2 of the License, or |
---|
| 26 | % (at your option) any later version. |
---|
| 27 | % |
---|
| 28 | % UVMAT is distributed in the hope that it will be useful, |
---|
| 29 | % but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 30 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 31 | % GNU General Public License (open UVMAT/COPYING.txt) for more details. |
---|
| 32 | %AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA |
---|
| 33 | |
---|
| 34 | function [Data,errormsg,result_conv]= civ_uvmat(Param,ncfile) |
---|
[246] | 35 | errormsg=''; |
---|
[242] | 36 | Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; |
---|
| 37 | Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes |
---|
| 38 | Data.Program='civ_uvmat'; |
---|
| 39 | Data.CivStage=0;%default |
---|
[270] | 40 | ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'}; %variables to read |
---|
[246] | 41 | ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'}; |
---|
[248] | 42 | mask=''; |
---|
| 43 | maskname='';%default |
---|
[251] | 44 | test_civx=0;%default |
---|
[270] | 45 | test_civ1=0;%default |
---|
| 46 | test_patch1=0;%default |
---|
[242] | 47 | |
---|
| 48 | %% Civ1 |
---|
| 49 | if isfield (Param,'Civ1') |
---|
[270] | 50 | test_civ1=1;% test for further use of civ1 results |
---|
[242] | 51 | par_civ1=Param.Civ1; |
---|
| 52 | image1=imread(par_civ1.filename_ima_a); |
---|
| 53 | image2=imread(par_civ1.filename_ima_b); |
---|
[309] | 54 | ibx2=ceil(par_civ1.Bx/2); |
---|
| 55 | iby2=ceil(par_civ1.By/2); |
---|
| 56 | isx2=ceil(par_civ1.Searchx/2); |
---|
| 57 | isy2=ceil(par_civ1.Searchy/2); |
---|
| 58 | shiftx=par_civ1.Shiftx; |
---|
| 59 | shifty=par_civ1.Shifty; |
---|
[246] | 60 | miniy=max(1+isy2+shifty,1+iby2); |
---|
[242] | 61 | minix=max(1+isx2-shiftx,1+ibx2); |
---|
[246] | 62 | maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2); |
---|
[242] | 63 | maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2); |
---|
[309] | 64 | if ~isfield(par_civ1,'PointCoord') |
---|
| 65 | [GridX,GridY]=meshgrid(minix:par_civ1.Dx:maxix,miniy:par_civ1.Dy:maxiy); |
---|
| 66 | par_civ1.PointCoord(:,1)=reshape(GridX,[],1); |
---|
| 67 | par_civ1.PointCoord(:,2)=reshape(GridY,[],1); |
---|
| 68 | end |
---|
| 69 | if par_civ1.CheckMask && isfield(par_civ1,'MaskName') && ~isempty(par_civ1.MaskName) |
---|
| 70 | maskname=par_civ1.MaskName; |
---|
[248] | 71 | mask=imread(maskname); |
---|
| 72 | end |
---|
[242] | 73 | % caluclate velocity data (y and v in indices, reverse to y component) |
---|
[309] | 74 | [xtable ytable utable vtable ctable F result_conv errormsg] = civ (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,par_civ1.PointCoord,par_civ1.Rho, mask); |
---|
[242] | 75 | list_param=(fieldnames(par_civ1))'; |
---|
| 76 | list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'}; |
---|
[309] | 77 | index_remove=zeros(size(list_param)); |
---|
| 78 | for name=list_remove %loop on the list of names |
---|
| 79 | index_remove=index_remove +strcmp(name{1},list_param);%index of the current name = name{1} |
---|
[242] | 80 | end |
---|
[309] | 81 | list_param(find(index_remove,1))=[]; |
---|
| 82 | Civ1_param=list_param;%initialisation |
---|
[242] | 83 | for ilist=1:length(list_param) |
---|
| 84 | Civ1_param{ilist}=['Civ1_' list_param{ilist}]; |
---|
[309] | 85 | Data.(['Civ1_' list_param{ilist}])=par_civ1.(list_param{ilist}); |
---|
[242] | 86 | end |
---|
| 87 | if isfield(Data,'Civ1_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile') |
---|
| 88 | Data.Civ1_gridname=''; |
---|
| 89 | end |
---|
| 90 | if isfield(Data,'Civ1_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile') |
---|
| 91 | Data.Civ1_maskname=''; |
---|
| 92 | end |
---|
| 93 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param {'Civ1_Time','Civ1_Dt'}]; |
---|
| 94 | Data.Civ1_Time=str2double(par_civ1.T0); |
---|
| 95 | Data.Civ1_Dt=str2double(par_civ1.Dt); |
---|
| 96 | Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};% cell array containing the names of the fields to record |
---|
| 97 | Data.VarDimName={'nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'}; |
---|
| 98 | Data.VarAttribute{1}.Role='coord_x'; |
---|
| 99 | Data.VarAttribute{2}.Role='coord_y'; |
---|
| 100 | Data.VarAttribute{3}.Role='vector_x'; |
---|
| 101 | Data.VarAttribute{4}.Role='vector_y'; |
---|
| 102 | Data.VarAttribute{5}.Role='warnflag'; |
---|
| 103 | Data.Civ1_X=reshape(xtable,[],1); |
---|
| 104 | Data.Civ1_Y=reshape(size(image1,1)-ytable+1,[],1); |
---|
| 105 | Data.Civ1_U=reshape(utable,[],1); |
---|
| 106 | Data.Civ1_V=reshape(-vtable,[],1); |
---|
| 107 | Data.Civ1_C=reshape(ctable,[],1); |
---|
| 108 | Data.Civ1_F=reshape(F,[],1); |
---|
| 109 | Data.CivStage=1; |
---|
| 110 | else |
---|
[270] | 111 | Data=nc2struct(ncfile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format |
---|
[251] | 112 | if ~isempty(Data.absolut_time_T0')%read civx file |
---|
[270] | 113 | test_civx=1;% test for old civx data format |
---|
[251] | 114 | [Data,vardetect,ichoice]=nc2struct(ncfile);%read the variables in the netcdf file |
---|
[246] | 115 | else |
---|
[251] | 116 | if isfield(Param,'Fix1') |
---|
| 117 | Data=nc2struct(ncfile,ListVarCiv1);%read civ1 data in the existing netcdf file |
---|
| 118 | else |
---|
| 119 | Data=nc2struct(ncfile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file |
---|
| 120 | end |
---|
[246] | 121 | end |
---|
| 122 | if isfield(Data,'Txt') |
---|
| 123 | msgbox_uvmat('ERROR',Data.Txt) |
---|
| 124 | return |
---|
| 125 | end |
---|
[242] | 126 | end |
---|
| 127 | |
---|
| 128 | %% Fix1 |
---|
[309] | 129 | % if Param.CheckCiv1==1 |
---|
| 130 | % Param.Civ1=Param.Civ1; |
---|
| 131 | % end |
---|
| 132 | % if Param.CheckFix1==1 |
---|
| 133 | % Param.Fix1=Param.Fix1; |
---|
| 134 | % fix1.WarnFlags=[]; |
---|
| 135 | % if get(handles.CheckFmin2,'Value') |
---|
| 136 | % fix1.WarnFlags=[fix1.WarnFlags -2]; |
---|
| 137 | % end |
---|
| 138 | % if get(handles.CheckF3,'Value') |
---|
| 139 | % fix1.WarnFlags=[fix1.WarnFlags 3]; |
---|
| 140 | % end |
---|
| 141 | % fix1.LowerBoundCorr=thresh_vecC1; |
---|
| 142 | % if get(handles.num_MinVel,'Value') |
---|
| 143 | % fix1.UppperBoundVel=thresh_vel1; |
---|
| 144 | % else |
---|
| 145 | % fix1.LowerBoundVel=thresh_vel1; |
---|
| 146 | % end |
---|
| 147 | % if get(handles.CheckMask,'Value') |
---|
| 148 | % fix1.MaskName=maskname; |
---|
| 149 | % end |
---|
| 150 | % Param.Fix1=fix1; |
---|
| 151 | % end |
---|
| 152 | % if Param.CheckPatch1==1 |
---|
| 153 | % if strcmp(compare,'stereo PIV') |
---|
| 154 | % filebase_A=filecell.filebase; |
---|
| 155 | % [pp,ff]=fileparts(filebase_A); |
---|
| 156 | % filebase_B=fullfile(pp,get(handles.RootName_1,'String')); |
---|
| 157 | % %TO CHECK: filecell.nc.civ1{ifile,j},filecell.ncA.civ1{ifile,j} have been switched according to Matias Duran |
---|
| 158 | % RUN_STLIN(filecell.nc.civ1{ifile,j},filecell.ncA.civ1{ifile,j},'civ1',filecell.st{ifile,j},... |
---|
| 159 | % str2num(nx_patch1),str2num(ny_patch1),str2num(thresh_patch1),[filebase_A '.xml'],[filebase_B '.xml']) |
---|
| 160 | % else |
---|
| 161 | % Param.Patch1.Rho=rho_patch1; |
---|
| 162 | % Param.Patch1.Threshold=thresh_patch1; |
---|
| 163 | % Param.Patch1.SubDomain=subdomain_patch1; |
---|
| 164 | % end |
---|
| 165 | % end |
---|
| 166 | % if Param.CheckCiv2==1 |
---|
| 167 | % Param.Civ2=Param.Civ2; |
---|
| 168 | % end |
---|
| 169 | % if Param.CheckFix2==1 |
---|
| 170 | % fix2.WarnFlags=[]; |
---|
| 171 | % if get(handles.CheckFmin2,'Value') |
---|
| 172 | % fix2.WarnFlags=[fix2.WarnFlags -2]; |
---|
| 173 | % end |
---|
| 174 | % if get(handles.CheckF4,'Value') |
---|
| 175 | % fix2.WarnFlags=[fix2.WarnFlags 4]; |
---|
| 176 | % end |
---|
| 177 | % if get(handles.CheckF3,'Value') |
---|
| 178 | % fix2.WarnFlags=[fix2.WarnFlags 3]; |
---|
| 179 | % end |
---|
| 180 | % fix2.LowerBoundCorr=thresh_vec2C; |
---|
| 181 | % if get(handles.num_MinVel,'Value') |
---|
| 182 | % fix2.UppperBoundVel=thresh_vel2; |
---|
| 183 | % else |
---|
| 184 | % fix2.LowerBoundVel=thresh_vel2; |
---|
| 185 | % end |
---|
| 186 | % if get(handles.CheckMask,'Value') |
---|
| 187 | % fix2.MaskName=maskname; |
---|
| 188 | % end |
---|
| 189 | % Param.Fix2=fix2; |
---|
| 190 | % end |
---|
| 191 | % if Param.CheckPatch2==1 |
---|
| 192 | % if strcmp(compare,'stereo PIV') |
---|
| 193 | % filebase_A=filecell.filebase; |
---|
| 194 | % [pp,ff]=fileparts(filebase_A); |
---|
| 195 | % filebase_B=fullfile(pp,get(handles.RootName_1,'String')); |
---|
| 196 | % RUN_STLIN(filecell.ncA.civ2{ifile,j},filecell.nc.civ2{ifile,j},'civ2',filecell.st{ifile,j},... |
---|
| 197 | % str2num(nx_patch2),str2num(ny_patch2),str2num(thresh_patch2),[filebase_A '.xml'],[filebase_B '.xml']) |
---|
| 198 | % else |
---|
| 199 | % Param.Patch2.Rho=rho_patch2; |
---|
| 200 | % Param.Patch2.Threshold=thresh_patch2; |
---|
| 201 | % Param.Patch2.SubDomain=subdomain_patch2; |
---|
| 202 | % end |
---|
| 203 | % end |
---|
| 204 | |
---|
[242] | 205 | if isfield (Param,'Fix1') |
---|
| 206 | ListFixParam=fieldnames(Param.Fix1); |
---|
| 207 | for ilist=1:length(ListFixParam) |
---|
| 208 | ParamName=ListFixParam{ilist}; |
---|
| 209 | ListName=['Fix1_' ParamName]; |
---|
| 210 | eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']) |
---|
| 211 | eval(['Data.' ListName '=Param.Fix1.' ParamName ';']) |
---|
| 212 | end |
---|
| 213 | % Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Fix1_WarnFlags','Fix1_TreshCorr','Fix1_TreshVel','Fix1_UpperBoundTest'}]; |
---|
| 214 | % Data.Fix1_WarnFlags=Param.Fix1.WarnFlags; |
---|
| 215 | % Data.Fix1_ThreshCorr=Param.Fix1.ThreshCorr; |
---|
| 216 | % Data.Fix1_ThreshVel=Param.Fix1.ThreshVel; |
---|
| 217 | % Data.Fix1_UpperBoundTest=Param.Fix1.UpperBoundTest; |
---|
[251] | 218 | |
---|
| 219 | if test_civx |
---|
| 220 | if ~isfield(Data,'fix') |
---|
| 221 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix']; |
---|
| 222 | Data.fix=1; |
---|
| 223 | Data.ListVarName=[Data.ListVarName {'vec_FixFlag'}]; |
---|
[252] | 224 | Data.VarDimName=[Data.VarDimName {'nb_vectors'}]; |
---|
[251] | 225 | end |
---|
[309] | 226 | Data.vec_FixFlag=fix(Param.Fix1,Data.vec_F,Data.vec_C,Data.vec_U,Data.vec_V,Data.vec_X,Data.vec_Y); |
---|
[251] | 227 | else |
---|
| 228 | Data.ListVarName=[Data.ListVarName {'Civ1_FF'}]; |
---|
| 229 | Data.VarDimName=[Data.VarDimName {'nbvec1'}]; |
---|
| 230 | nbvar=length(Data.ListVarName); |
---|
| 231 | Data.VarAttribute{nbvar}.Role='errorflag'; |
---|
[309] | 232 | Data.Civ1_FF=fix(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V); |
---|
[251] | 233 | Data.CivStage=2; |
---|
| 234 | end |
---|
[242] | 235 | end |
---|
| 236 | %% Patch1 |
---|
| 237 | if isfield (Param,'Patch1') |
---|
[270] | 238 | test_patch1=1; |
---|
[242] | 239 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}]; |
---|
| 240 | Data.Patch1_Rho=str2double(Param.Patch1.Rho); |
---|
| 241 | Data.Patch1_Threshold=str2double(Param.Patch1.Threshold); |
---|
| 242 | Data.Patch1_SubDomain=str2double(Param.Patch1.SubDomain); |
---|
[246] | 243 | Data.ListVarName=[Data.ListVarName {'Civ1_U_Diff','Civ1_V_Diff','Civ1_X_SubRange','Civ1_Y_SubRange','Civ1_NbSites','Civ1_X_tps','Civ1_Y_tps','Civ1_U_tps','Civ1_V_tps'}]; |
---|
| 244 | Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbSubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',... |
---|
| 245 | {'NbVec1Sub','NbSubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}]; |
---|
[242] | 246 | nbvar=length(Data.ListVarName); |
---|
| 247 | Data.VarAttribute{nbvar-1}.Role='vector_x'; |
---|
| 248 | Data.VarAttribute{nbvar}.Role='vector_y'; |
---|
[246] | 249 | Data.Civ1_U_Diff=zeros(size(Data.Civ1_X)); |
---|
| 250 | Data.Civ1_V_Diff=zeros(size(Data.Civ1_X)); |
---|
[242] | 251 | if isfield(Data,'Civ1_FF') |
---|
| 252 | ind_good=find(Data.Civ1_FF==0); |
---|
| 253 | else |
---|
| 254 | ind_good=1:numel(Data.Civ1_X); |
---|
| 255 | end |
---|
[246] | 256 | [Data.Civ1_X_SubRange,Data.Civ1_Y_SubRange,Data.Civ1_NbSites,FFres,Ures, Vres,Data.Civ1_X_tps,Data.Civ1_Y_tps,Data.Civ1_U_tps,Data.Civ1_V_tps]=... |
---|
[309] | 257 | patch(Data.Civ1_X(ind_good)',Data.Civ1_Y(ind_good)',Data.Civ1_U(ind_good)',Data.Civ1_V(ind_good)',Data.Patch1_Rho,Data.Patch1_Threshold,Data.Patch1_SubDomain); |
---|
[246] | 258 | Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)-Ures; |
---|
| 259 | Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)-Vres; |
---|
| 260 | Data.Civ1_FF(ind_good)=FFres; |
---|
| 261 | Data.CivStage=3; |
---|
[242] | 262 | end |
---|
[248] | 263 | |
---|
| 264 | %% Civ2 |
---|
| 265 | if isfield (Param,'Civ2') |
---|
[252] | 266 | par_civ2=Param.Civ2; |
---|
[270] | 267 | if ~test_civ1 || ~strcmp(par_civ1.filename_ima_a,par_civ2.filename_ima_a) |
---|
| 268 | image1=imread(par_civ2.filename_ima_a);%read first image if not already done for civ1 |
---|
| 269 | end |
---|
| 270 | if ~test_civ1|| ~strcmp(par_civ1.filename_ima_b,par_civ2.filename_ima_b) |
---|
| 271 | image2=imread(par_civ2.filename_ima_b);%read second image if not already done for civ1 |
---|
| 272 | end |
---|
[309] | 273 | % stepx=str2double(par_civ2.dx); |
---|
| 274 | % stepy=str2double(par_civ2.dy); |
---|
[270] | 275 | ibx2=ceil(str2double(par_civ2.ibx)/2); |
---|
| 276 | iby2=ceil(str2double(par_civ2.iby)/2); |
---|
| 277 | isx2=ibx2+2; |
---|
| 278 | isy2=iby2+2; |
---|
| 279 | %get the previous guess for displacement |
---|
| 280 | if ~test_patch1 |
---|
| 281 | [Guess,VelType]=read_civdata(ObjectName,InputField,ParamIn.VelType);%TO DEVELOP |
---|
| 282 | else |
---|
| 283 | Data.Civ1_U_Diff=zeros(size(Data.Civ1_X)); |
---|
| 284 | Data.Civ1_V_Diff=zeros(size(Data.Civ1_X)); |
---|
| 285 | if isfield(Data,'Civ1_FF') |
---|
| 286 | ind_good=find(Data.Civ1_FF==0); |
---|
| 287 | else |
---|
| 288 | ind_good=1:numel(Data.Civ1_X); |
---|
| 289 | end |
---|
| 290 | [Data.Civ1_X_SubRange,Data.Civ1_Y_SubRange,Data.Civ1_NbSites,FFres,Ures, Vres,Data.Civ1_X_tps,Data.Civ1_Y_tps,Data.Civ1_U_tps,Data.Civ1_V_tps]=... |
---|
[309] | 291 | patch(Data.Civ1_X(ind_good)',Data.Civ1_Y(ind_good)',Data.Civ1_U(ind_good)',Data.Civ1_V(ind_good)',Data.Patch1_Rho,Data.Patch1_Threshold,Data.Patch1_SubDomain); |
---|
[270] | 292 | end |
---|
[248] | 293 | % shiftx=str2num(par_civ1.shiftx); |
---|
| 294 | % shifty=str2num(par_civ1.shifty); |
---|
| 295 | % TO GET shift from par_civ2.filename_nc1 |
---|
[270] | 296 | % shiftx=velocity interpolated at position |
---|
[248] | 297 | miniy=max(1+isy2+shifty,1+iby2); |
---|
| 298 | minix=max(1+isx2-shiftx,1+ibx2); |
---|
| 299 | maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2); |
---|
| 300 | maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2); |
---|
[309] | 301 | [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy); |
---|
[248] | 302 | PointCoord(:,1)=reshape(GridX,[],1); |
---|
| 303 | PointCoord(:,2)=reshape(GridY,[],1); |
---|
| 304 | if ~isempty(par_civ2.maskname)&& ~strcmp(maskname,par_civ2.maskname)% mask exist, not already read in civ1 |
---|
| 305 | mask=imread(par_civ2.maskname); |
---|
| 306 | end |
---|
| 307 | % caluclate velocity data (y and v in indices, reverse to y component) |
---|
[309] | 308 | [xtable ytable utable vtable ctable F] = civ (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho),mask); |
---|
[248] | 309 | list_param=(fieldnames(par_civ1))'; |
---|
| 310 | list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'}; |
---|
| 311 | index=zeros(size(list_param)); |
---|
| 312 | for ilist=1:length(list_remove) |
---|
| 313 | index=strcmp(list_remove{ilist},list_param); |
---|
| 314 | if ~isempty(find(index,1)) |
---|
| 315 | list_param(index)=[]; |
---|
| 316 | end |
---|
| 317 | end |
---|
| 318 | for ilist=1:length(list_param) |
---|
| 319 | Civ1_param{ilist}=['Civ1_' list_param{ilist}]; |
---|
| 320 | eval(['Data.Civ1_' list_param{ilist} '=Param.Civ1.' list_param{ilist} ';']) |
---|
| 321 | end |
---|
| 322 | if isfield(Data,'Civ1_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile') |
---|
| 323 | Data.Civ1_gridname=''; |
---|
| 324 | end |
---|
| 325 | if isfield(Data,'Civ1_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile') |
---|
| 326 | Data.Civ1_maskname=''; |
---|
| 327 | end |
---|
| 328 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param {'Civ1_Time','Civ1_Dt'}]; |
---|
| 329 | Data.Civ1_Time=str2double(par_civ1.T0); |
---|
| 330 | Data.Civ1_Dt=str2double(par_civ1.Dt); |
---|
| 331 | Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};% cell array containing the names of the fields to record |
---|
| 332 | Data.VarDimName={'nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'}; |
---|
| 333 | Data.VarAttribute{1}.Role='coord_x'; |
---|
| 334 | Data.VarAttribute{2}.Role='coord_y'; |
---|
| 335 | Data.VarAttribute{3}.Role='vector_x'; |
---|
| 336 | Data.VarAttribute{4}.Role='vector_y'; |
---|
| 337 | Data.VarAttribute{5}.Role='warnflag'; |
---|
| 338 | Data.Civ1_X=reshape(xtable,[],1); |
---|
| 339 | Data.Civ1_Y=reshape(size(image1,1)-ytable+1,[],1); |
---|
| 340 | Data.Civ1_U=reshape(utable,[],1); |
---|
| 341 | Data.Civ1_V=reshape(-vtable,[],1); |
---|
| 342 | Data.Civ1_C=reshape(ctable,[],1); |
---|
| 343 | Data.Civ1_F=reshape(F,[],1); |
---|
| 344 | Data.CivStage=Data.CivStage+1; |
---|
| 345 | end |
---|
[252] | 346 | |
---|
| 347 | %% Fix2 |
---|
| 348 | if isfield (Param,'Fix2') |
---|
| 349 | ListFixParam=fieldnames(Param.Fix2); |
---|
| 350 | for ilist=1:length(ListFixParam) |
---|
| 351 | ParamName=ListFixParam{ilist}; |
---|
| 352 | ListName=['Fix1_' ParamName]; |
---|
| 353 | eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']) |
---|
| 354 | eval(['Data.' ListName '=Param.Fix2.' ParamName ';']) |
---|
| 355 | end |
---|
| 356 | if test_civx |
---|
| 357 | if ~isfield(Data,'fix2') |
---|
| 358 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2']; |
---|
| 359 | Data.fix2=1; |
---|
| 360 | Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}]; |
---|
| 361 | Data.VarDimName=[Data.VarDimName {'nb_vectors2'}]; |
---|
| 362 | end |
---|
[309] | 363 | Data.vec_FixFlag=fix(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y); |
---|
[252] | 364 | else |
---|
| 365 | Data.ListVarName=[Data.ListVarName {'Civ2_FF'}]; |
---|
| 366 | Data.VarDimName=[Data.VarDimName {'nbvec2'}]; |
---|
| 367 | nbvar=length(Data.ListVarName); |
---|
| 368 | Data.VarAttribute{nbvar}.Role='errorflag'; |
---|
[309] | 369 | Data.Civ2_FF=fix(Param.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V); |
---|
[252] | 370 | Data.CivStage=5; |
---|
| 371 | end |
---|
[270] | 372 | |
---|
[252] | 373 | end |
---|
| 374 | |
---|
| 375 | %% Patch2 |
---|
| 376 | if isfield (Param,'Patch2') |
---|
| 377 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch2_Rho','Patch2_Threshold','Patch2_SubDomain'}]; |
---|
| 378 | Data.Patch2_Rho=str2double(Param.Patch2.Rho); |
---|
| 379 | Data.Patch2_Threshold=str2double(Param.Patch2.Threshold); |
---|
| 380 | Data.Patch2_SubDomain=str2double(Param.Patch2.SubDomain); |
---|
| 381 | Data.ListVarName=[Data.ListVarName {'Civ2_U_Diff','Civ2_V_Diff','Civ2_X_SubRange','Civ2_Y_SubRange','Civ2_NbSites','Civ2_X_tps','Civ2_Y_tps','Civ2_U_tps','Civ2_V_tps'}]; |
---|
| 382 | Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2',{'NbSubDomain2','Two'},{'NbSubDomain2','Two'},'NbSubDomain2',... |
---|
| 383 | {'NbVec2Sub','NbSubDomain2'},{'NbVec2Sub','NbSubDomain2'},{'Nbtps2','NbSubDomain2'},{'Nbtps2','NbSubDomain2'}}]; |
---|
| 384 | nbvar=length(Data.ListVarName); |
---|
| 385 | Data.VarAttribute{nbvar-1}.Role='vector_x'; |
---|
| 386 | Data.VarAttribute{nbvar}.Role='vector_y'; |
---|
| 387 | Data.Civ2_U_Diff=zeros(size(Data.Civ2_X)); |
---|
| 388 | Data.Civ2_V_Diff=zeros(size(Data.Civ2_X)); |
---|
| 389 | if isfield(Data,'Civ2_FF') |
---|
| 390 | ind_good=find(Data.Civ2_FF==0); |
---|
| 391 | else |
---|
| 392 | ind_good=1:numel(Data.Civ2_X); |
---|
| 393 | end |
---|
| 394 | [Data.Civ2_X_SubRange,Data.Civ2_Y_SubRange,Data.Civ2_NbSites,FFres,Ures, Vres,Data.Civ2_X_tps,Data.Civ2_Y_tps,Data.Civ2_U_tps,Data.Civ2_V_tps]=... |
---|
[309] | 395 | patch(Data.Civ2_X(ind_good)',Data.Civ2_Y(ind_good)',Data.Civ2_U(ind_good)',Data.Civ2_V(ind_good)',Data.Patch2_Rho,Data.Patch2_Threshold,Data.Patch2_SubDomain); |
---|
[252] | 396 | Data.Civ2_U_Diff(ind_good)=Data.Civ2_U(ind_good)-Ures; |
---|
| 397 | Data.Civ2_V_Diff(ind_good)=Data.Civ2_V(ind_good)-Vres; |
---|
| 398 | Data.Civ2_FF(ind_good)=FFres; |
---|
[270] | 399 | Data.CivStage=6; |
---|
[252] | 400 | end |
---|
| 401 | |
---|
[309] | 402 | %% write result in a netcdf file if requested |
---|
[246] | 403 | if exist('ncfile','var') |
---|
[309] | 404 | errormsg=struct2nc(ncfile,Data); |
---|
[246] | 405 | end |
---|
[242] | 406 | |
---|
| 407 | |
---|
| 408 | %'RUN_FIX': function for fixing velocity fields: |
---|
| 409 | %----------------------------------------------- |
---|
| 410 | % RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref) |
---|
| 411 | % |
---|
| 412 | %filename: name of the netcdf file (used as input and output) |
---|
| 413 | %field: structure specifying the names of the fields to fix (depending on civ1 or civ2) |
---|
| 414 | %.vel_type='civ1' or 'civ2'; |
---|
| 415 | %.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1); |
---|
| 416 | %.fixflag=name of fix flag variable ('vec_FixFlag' for civ1) |
---|
| 417 | %flagindex: flag specifying which values of vec_f are removed: |
---|
| 418 | % if flagindex(1)=1: vec_f=-2 vectors are removed |
---|
| 419 | % if flagindex(2)=1: vec_f=3 vectors are removed |
---|
| 420 | % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2) |
---|
| 421 | %iter=1 for civ1 fields and iter=2 for civ2 fields |
---|
| 422 | %thresh_vecC: threshold in the image correlation vec_C |
---|
| 423 | %flag_mask: =1 mask used to remove vectors (0 else) |
---|
| 424 | %maskname: name of the mask image file for fix |
---|
| 425 | %thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists |
---|
| 426 | %inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold |
---|
| 427 | %fileref: .nc file name for a reference velocity (='': refrence 0 used) |
---|
| 428 | %fieldref: 'civ1','filter1'...feld used in fileref |
---|
| 429 | |
---|
[309] | 430 | function FF=fix(Param,F,C,U,V,X,Y) |
---|
[242] | 431 | FF=zeros(size(F));%default |
---|
[309] | 432 | Param |
---|
[242] | 433 | |
---|
| 434 | %criterium on warn flags |
---|
[309] | 435 | FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'}; |
---|
| 436 | FlagVal=[-2 2 3 4]; |
---|
| 437 | for iflag=1:numel(FlagVal) |
---|
| 438 | if Param.(CheckFlag(iflag)) |
---|
| 439 | FF=(FF==1| F==FlagVal(iflag)); |
---|
| 440 | % if isfield (Param,'WarnFlags') |
---|
| 441 | % for iflag=1:numel(Param.WarnFlags) |
---|
| 442 | % FF=(FF==1| F==Param.WarnFlags(iflag)); |
---|
| 443 | % end |
---|
| 444 | % end |
---|
[242] | 445 | end |
---|
| 446 | end |
---|
| 447 | %criterium on correlation values |
---|
[309] | 448 | if isfield (Param,'MinCorr') |
---|
| 449 | FF=FF==1 | C<Param.MinCorr; |
---|
[242] | 450 | end |
---|
[309] | 451 | return |
---|
[242] | 452 | |
---|
[309] | 453 | |
---|
[242] | 454 | if isfield (Param,'LowerBoundVel')&& ~isequal(Param.LowerBoundVel,0) |
---|
| 455 | thresh=Param.LowerBoundVel*Param.LowerBoundVel; |
---|
| 456 | FF=FF==1 | (U.*U+V.*V)<thresh; |
---|
| 457 | end |
---|
| 458 | if isfield (Param,'UpperBoundVel')&& ~isequal(Param.UpperBoundVel,0) |
---|
| 459 | thresh=Param.UpperBoundVel*Param.UpperBoundVel; |
---|
| 460 | FF=FF==1 | (U.*U+V.*V)>thresh; |
---|
| 461 | end |
---|
[251] | 462 | if isfield(Param,'MaskName') |
---|
| 463 | M=imread(Param.MaskName); |
---|
| 464 | nxy=size(M); |
---|
| 465 | M=reshape(M,1,[]); |
---|
| 466 | rangx0=[0.5 nxy(2)-0.5]; |
---|
| 467 | rangy0=[0.5 nxy(1)-0.5]; |
---|
| 468 | vec_x1=X-U/2;%beginning points |
---|
| 469 | vec_x2=X+U/2;%end points of vectors |
---|
| 470 | vec_y1=Y-V/2;%beginning points |
---|
| 471 | vec_y2=Y+V/2;%end points of vectors |
---|
| 472 | indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x1 |
---|
| 473 | indy=1+round((nxy(1)-1)*(vec_y1-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y1 |
---|
| 474 | test_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside |
---|
| 475 | indx=indx.*test_in+(1-test_in); %replace indx by 1 out of the mask range |
---|
| 476 | indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range |
---|
| 477 | ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector |
---|
| 478 | Mvalues=M(ICOMB); |
---|
| 479 | flag7b=((20 < Mvalues) & (Mvalues < 200))| ~test_in'; |
---|
| 480 | indx=1+round((nxy(2)-1)*(vec_x2-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x2 |
---|
| 481 | indy=1+round((nxy(1)-1)*(vec_y2-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y2 |
---|
| 482 | test_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside |
---|
| 483 | indx=indx.*test_in+(1-test_in); %replace indx by 1 out of the mask range |
---|
| 484 | indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range |
---|
| 485 | ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector |
---|
| 486 | Mvalues=M(ICOMB); |
---|
| 487 | flag7e=((Mvalues > 20) & (Mvalues < 200))| ~test_in'; |
---|
| 488 | FF=FF==1 |(flag7b|flag7e)'; |
---|
| 489 | end |
---|
| 490 | % flag7=0; |
---|
| 491 | % end |
---|
| 492 | |
---|
| 493 | |
---|
[246] | 494 | FF=double(FF); |
---|
[242] | 495 | % |
---|
| 496 | % % criterium on velocity values |
---|
| 497 | % delta_u=Field.U;%default without ref file |
---|
| 498 | % delta_v=Field.V; |
---|
| 499 | % if exist('fileref','var') && ~isempty(fileref) |
---|
| 500 | % if ~exist(fileref,'file') |
---|
| 501 | % error='reference file not found in RUN_FIX.m'; |
---|
| 502 | % display(error); |
---|
| 503 | % return |
---|
| 504 | % end |
---|
| 505 | % FieldRef=read_civxdata(fileref,[],fieldref); |
---|
| 506 | % if isfield(FieldRef,'FF') |
---|
| 507 | % index_true=find(FieldRef.FF==0); |
---|
| 508 | % FieldRef.X=FieldRef.X(index_true); |
---|
| 509 | % FieldRef.Y=FieldRef.Y(index_true); |
---|
| 510 | % FieldRef.U=FieldRef.U(index_true); |
---|
| 511 | % FieldRef.V=FieldRef.V(index_true); |
---|
| 512 | % end |
---|
| 513 | % if ~isfield(FieldRef,'X') || ~isfield(FieldRef,'Y') || ~isfield(FieldRef,'U') || ~isfield(FieldRef,'V') |
---|
| 514 | % error='reference file is not a velocity field in RUN_FIX.m '; %bad input file |
---|
| 515 | % return |
---|
| 516 | % end |
---|
| 517 | % if length(FieldRef.X)<=1 |
---|
| 518 | % errordlg('reference field with one vector or less in RUN_FIX.m') |
---|
| 519 | % return |
---|
| 520 | % end |
---|
| 521 | % vec_U_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.U,Field.X,Field.Y); %interpolate vectors in the ref field |
---|
| 522 | % vec_V_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.V,Field.X,Field.Y); %interpolate vectors in the ref field to the positions of the main field |
---|
| 523 | % delta_u=Field.U-vec_U_ref;%take the difference with the interpolated ref field |
---|
| 524 | % delta_v=Field.V-vec_V_ref; |
---|
| 525 | % end |
---|
| 526 | % thresh_vel_x=thresh_vel; |
---|
| 527 | % thresh_vel_y=thresh_vel; |
---|
| 528 | % if isequal(inf_sup,1) |
---|
| 529 | % flag5=abs(delta_u)<thresh_vel_x & abs(delta_v)<thresh_vel_y &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1); |
---|
| 530 | % elseif isequal(inf_sup,2) |
---|
| 531 | % flag5=(abs(delta_u)>thresh_vel_x | abs(delta_v)>thresh_vel_y) &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1); |
---|
| 532 | % end |
---|
| 533 | % |
---|
| 534 | % % flag7 introduce a grey mask, matrix M |
---|
[251] | 535 | |
---|
[242] | 536 | % flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7; |
---|
| 537 | % fixflag_unit=Field.FF-10*floor(Field.FF/10); %unity term of fix_flag |
---|
| 538 | |
---|
| 539 | |
---|
| 540 | |
---|
| 541 | %------------------------------------------------------------------------ |
---|
| 542 | % patch function |
---|
[309] | 543 | function [SubRangx,SubRangy,nbpoints,FF,U_smooth,V_smooth,X_tps,Y_tps,U_tps,V_tps] =patch(X,Y,U,V,Rho,Threshold,SubDomain) |
---|
[242] | 544 | %subdomain decomposition |
---|
| 545 | warning off |
---|
[243] | 546 | U=reshape(U,[],1); |
---|
| 547 | V=reshape(V,[],1); |
---|
| 548 | X=reshape(X,[],1); |
---|
| 549 | Y=reshape(Y,[],1); |
---|
[242] | 550 | nbvec=numel(X); |
---|
[246] | 551 | NbSubDomain=ceil(nbvec/SubDomain); |
---|
| 552 | MinX=min(X); |
---|
| 553 | MinY=min(Y); |
---|
| 554 | MaxX=max(X); |
---|
| 555 | MaxY=max(Y); |
---|
[242] | 556 | RangX=MaxX-MinX; |
---|
| 557 | RangY=MaxY-MinY; |
---|
[246] | 558 | AspectRatio=RangY/RangX; |
---|
| 559 | NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1); |
---|
| 560 | NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1); |
---|
| 561 | NbSubDomain=NbSubDomainX*NbSubDomainY; |
---|
[242] | 562 | SizX=RangX/NbSubDomainX;%width of subdomains |
---|
| 563 | SizY=RangY/NbSubDomainY;%height of subdomains |
---|
| 564 | CentreX=linspace(MinX+SizX/2,MaxX-SizX/2,NbSubDomainX); |
---|
| 565 | CentreY=linspace(MinY+SizY/2,MaxY-SizY/2,NbSubDomainY); |
---|
[246] | 566 | [CentreX,CentreY]=meshgrid(CentreX,CentreY); |
---|
| 567 | CentreY=reshape(CentreY,1,[]); |
---|
| 568 | CentreX=reshape(CentreX,1,[]); |
---|
| 569 | rho=SizX*SizY*Rho/1000000;%optimum rho increase as the area of the subdomain (division by 10^6 to reach good values with the default GUI input) |
---|
| 570 | U_tps_sub=zeros(length(X),NbSubDomain);%default spline |
---|
| 571 | V_tps_sub=zeros(length(X),NbSubDomain);%default spline |
---|
| 572 | U_smooth=zeros(length(X),1); |
---|
| 573 | V_smooth=zeros(length(X),1); |
---|
| 574 | |
---|
| 575 | nb_select=zeros(length(X),1); |
---|
[243] | 576 | FF=zeros(length(X),1); |
---|
[246] | 577 | test_empty=zeros(1,NbSubDomain); |
---|
| 578 | for isub=1:NbSubDomain |
---|
| 579 | SubRangx(isub,:)=[CentreX(isub)-SizX/2 CentreX(isub)+SizX/2]; |
---|
| 580 | SubRangy(isub,:)=[CentreY(isub)-SizY/2 CentreY(isub)+SizY/2]; |
---|
| 581 | ind_sel_previous=[]; |
---|
| 582 | ind_sel=0; |
---|
| 583 | while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most |
---|
| 584 | ind_sel_previous=ind_sel; |
---|
| 585 | ind_sel=find(X>SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2)); |
---|
| 586 | % if no vector in the subdomain, skip the subdomain |
---|
| 587 | if isempty(ind_sel) |
---|
| 588 | test_empty(isub)=1; |
---|
| 589 | U_tps(1,isub)=0;%define U_tps and V_tps by default |
---|
| 590 | V_tps(1,isub)=0; |
---|
| 591 | break |
---|
| 592 | % if too few selected vectors, increase the subrange for next iteration |
---|
| 593 | elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); |
---|
| 594 | SubRangx(isub,1)=SubRangx(isub,1)-SizX/4; |
---|
| 595 | SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; |
---|
| 596 | SubRangy(isub,1)=SubRangy(isub,1)-SizY/4; |
---|
| 597 | SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; |
---|
| 598 | else |
---|
| 599 | [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),U(ind_sel),rho); |
---|
| 600 | [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),V(ind_sel),rho); |
---|
| 601 | UDiff=U_smooth_sub-U(ind_sel); |
---|
| 602 | VDiff=V_smooth_sub-V(ind_sel); |
---|
| 603 | NormDiff=UDiff.*UDiff+VDiff.*VDiff; |
---|
| 604 | FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination |
---|
| 605 | ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors |
---|
| 606 | % no value exceeds threshold, the result is recorded |
---|
| 607 | if isequal(numel(ind_ind_sel),numel(ind_sel)) |
---|
| 608 | U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub; |
---|
| 609 | V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub; |
---|
| 610 | nbpoints(isub)=numel(ind_sel); |
---|
| 611 | X_tps(1:nbpoints(isub),isub)=X(ind_sel); |
---|
| 612 | Y_tps(1:nbpoints(isub),isub)=Y(ind_sel); |
---|
| 613 | U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; |
---|
| 614 | V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; |
---|
| 615 | nb_select(ind_sel)=nb_select(ind_sel)+1; |
---|
| 616 | display('good') |
---|
| 617 | break |
---|
| 618 | % too few selected vectors, increase the subrange for next iteration |
---|
| 619 | elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); |
---|
| 620 | SubRangx(isub,1)=SubRangx(isub,1)-SizX/4; |
---|
| 621 | SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; |
---|
| 622 | SubRangy(isub,1)=SubRangy(isub,1)-SizY/4; |
---|
| 623 | SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; |
---|
| 624 | % display('fewsmooth') |
---|
| 625 | % interpolation-smoothing is done again with the selected vectors |
---|
[242] | 626 | else |
---|
[246] | 627 | [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),U(ind_sel(ind_ind_sel)),rho); |
---|
| 628 | [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel(ind_ind_sel)),Y(ind_sel(ind_ind_sel)),V(ind_sel(ind_ind_sel)),rho); |
---|
| 629 | U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub; |
---|
| 630 | V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub; |
---|
| 631 | nbpoints(isub)=numel(ind_ind_sel); |
---|
| 632 | X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel)); |
---|
| 633 | Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel)); |
---|
| 634 | U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; |
---|
| 635 | V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; |
---|
| 636 | nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1; |
---|
| 637 | display('good2') |
---|
| 638 | break |
---|
| 639 | end |
---|
| 640 | end |
---|
[242] | 641 | end |
---|
| 642 | end |
---|
[246] | 643 | ind_empty=find(test_empty); |
---|
| 644 | %remove empty subdomains |
---|
| 645 | if ~isempty(ind_empty) |
---|
| 646 | SubRangx(ind_empty,:)=[]; |
---|
| 647 | SubRangy(ind_empty,:)=[]; |
---|
| 648 | X_tps(:,ind_empty)=[]; |
---|
| 649 | Y_tps(:,ind_empty)=[]; |
---|
| 650 | U_tps(:,ind_empty)=[]; |
---|
| 651 | V_tps(:,ind_empty)=[]; |
---|
| 652 | end |
---|
| 653 | nb_select(nb_select==0)=1;%ones(size(find(nb_select==0))); |
---|
| 654 | U_smooth=U_smooth./nb_select; |
---|
| 655 | V_smooth=V_smooth./nb_select; |
---|
[242] | 656 | |
---|
| 657 | |
---|
[309] | 658 | % 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/ |
---|
| 659 | %-------------------------------------------------------------------------- |
---|
| 660 | % function [xtable ytable utable vtable typevector] = civ (image1,image2,ibx,iby step, subpixfinder, mask, roi) |
---|
| 661 | % |
---|
| 662 | % OUTPUT: |
---|
| 663 | % xtable: set of x coordinates |
---|
| 664 | % ytable: set of y coordiantes |
---|
| 665 | % utable: set of u displacements (along x) |
---|
| 666 | % vtable: set of v displacements (along y) |
---|
| 667 | % ctable: max image correlation for each vector |
---|
| 668 | % typevector: set of flags, =1 for good, =0 for NaN vectors |
---|
| 669 | % |
---|
| 670 | %INPUT: |
---|
| 671 | % image1:first image (matrix) |
---|
| 672 | % image2: second image (matrix) |
---|
| 673 | % ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1) |
---|
| 674 | % isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1) |
---|
| 675 | % shiftx, shifty: shift of the search box (in pixel index, yshift reversed) |
---|
| 676 | % step: mesh of the measurement points (in px) |
---|
| 677 | % subpixfinder=1 or 2 controls the curve fitting of the image correlation |
---|
| 678 | % mask: =[] for no mask |
---|
| 679 | % roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[]; |
---|
| 680 | function [xtable ytable utable vtable ctable F result_conv errormsg] = civ (image1,image2,ibx2,iby2,isx2,isy2,shiftx,shifty, GridIndices, subpixfinder,mask) |
---|
| 681 | %this funtion performs the DCC PIV analysis. Recent window-deformation |
---|
| 682 | %methods perform better and will maybe be implemented in the future. |
---|
| 683 | nbvec=size(GridIndices,1); |
---|
| 684 | xtable=zeros(nbvec,1); |
---|
| 685 | ytable=xtable; |
---|
| 686 | utable=xtable; |
---|
| 687 | vtable=xtable; |
---|
| 688 | ctable=xtable; |
---|
| 689 | F=xtable; |
---|
| 690 | result_conv=[]; |
---|
| 691 | errormsg=''; |
---|
| 692 | %warning off %MATLAB:log:logOfZero |
---|
| 693 | [npy_ima npx_ima]=size(image1); |
---|
| 694 | if ~isequal(size(image2),[npy_ima npx_ima]) |
---|
| 695 | errormsg='image pair with unequal size'; |
---|
| 696 | return |
---|
| 697 | end |
---|
[242] | 698 | |
---|
[309] | 699 | %% mask |
---|
| 700 | testmask=0; |
---|
| 701 | image1=double(image1); |
---|
| 702 | image2=double(image2); |
---|
| 703 | if exist('mask','var') && ~isempty(mask) |
---|
| 704 | testmask=1; |
---|
| 705 | if ~isequal(size(mask),[npy_ima npx_ima]) |
---|
| 706 | errormsg='mask must be an image with the same size as the images'; |
---|
| 707 | return |
---|
| 708 | end |
---|
| 709 | % Convention for mask |
---|
| 710 | % mask >200 : velocity calculated |
---|
| 711 | % 200 >=mask>150;velocity not calculated, interpolation allowed (bad spots) |
---|
| 712 | % 150>=mask >100: velocity not calculated, nor interpolated |
---|
| 713 | % 100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries) |
---|
| 714 | % 20>=mask: velocity=0 |
---|
| 715 | test_noflux=(mask<100) ; |
---|
| 716 | test_undefined=(mask<200 & mask>=100 ); |
---|
| 717 | image1(test_undefined)=min(min(image1));% put image to zero in the undefined area |
---|
| 718 | image2(test_undefined)=min(min(image2));% put image to zero in the undefined area |
---|
| 719 | end |
---|
[246] | 720 | |
---|
[309] | 721 | %% calculate correlations: MAINLOOP on velocity vectors |
---|
| 722 | corrmax=0; |
---|
| 723 | sum_square=1;% default |
---|
| 724 | for ivec=1:nbvec |
---|
| 725 | iref=GridIndices(ivec,1); |
---|
| 726 | jref=GridIndices(ivec,2); |
---|
| 727 | testmask_ij=0; |
---|
| 728 | test0=0; |
---|
| 729 | if testmask |
---|
| 730 | if mask(jref,iref)<=20 |
---|
| 731 | vector=[0 0]; |
---|
| 732 | test0=1; |
---|
| 733 | else |
---|
| 734 | mask_crop1=mask(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2); |
---|
| 735 | mask_crop2=mask(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2); |
---|
| 736 | if ~isempty(find(mask_crop1<=200 & mask_crop1>100,1)) || ~isempty(find(mask_crop2<=200 & mask_crop2>100,1)); |
---|
| 737 | testmask_ij=1; |
---|
| 738 | end |
---|
| 739 | end |
---|
| 740 | end |
---|
| 741 | if ~test0 |
---|
| 742 | image1_crop=image1(jref-iby2:jref+iby2,iref-ibx2:iref+ibx2);%extract a subimage (correlation box) from images 1 |
---|
| 743 | image2_crop=image2(jref+shifty-isy2:jref+shifty+isy2,iref+shiftx-isx2:iref+shiftx+isx2);%extract a larger subimage (search box) from image 2 |
---|
| 744 | image1_crop=image1_crop-mean(mean(image1_crop));%substract the mean |
---|
| 745 | image2_crop=image2_crop-mean(mean(image2_crop)); |
---|
| 746 | %reference: Oliver Pust, PIV: Direct Cross-Correlation |
---|
| 747 | result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid'); |
---|
| 748 | corrmax= max(max(result_conv)); |
---|
| 749 | result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255 |
---|
| 750 | %Find the correlation max, at 255 |
---|
| 751 | [y,x] = find(result_conv==255,1); |
---|
| 752 | if ~isempty(y) && ~isempty(x) |
---|
| 753 | try |
---|
| 754 | if subpixfinder==1 |
---|
| 755 | [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y); |
---|
| 756 | elseif subpixfinder==2 |
---|
| 757 | [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y); |
---|
| 758 | end |
---|
| 759 | sum_square=sum(sum(image1_crop.*image1_crop)); |
---|
| 760 | ctable(ivec)=corrmax/sum_square;% correlation value |
---|
| 761 | % if vector(1)>shiftx+isx2-ibx2+subpixfinder || vector(2)>shifty+isy2-iby2+subpixfinder |
---|
| 762 | % F(ivec)=-2;%vector reaches the border of the search zone |
---|
| 763 | % end |
---|
| 764 | catch ME |
---|
| 765 | vector=[0 0]; %if something goes wrong with cross correlation..... |
---|
| 766 | F(ivec)=3; |
---|
| 767 | end |
---|
| 768 | else |
---|
| 769 | vector=[0 0]; %if something goes wrong with cross correlation..... |
---|
| 770 | F(ivec)=3; |
---|
| 771 | end |
---|
| 772 | if testmask_ij |
---|
| 773 | F(ivec)=3; |
---|
| 774 | end |
---|
| 775 | end |
---|
| 776 | |
---|
| 777 | %Create the vector matrix x, y, u, v |
---|
| 778 | xtable(ivec)=iref+vector(1)/2;% convec flow (velocity taken at the point middle from imgae1 and 2) |
---|
| 779 | ytable(ivec)=jref+vector(2)/2; |
---|
| 780 | utable(ivec)=vector(1)+shiftx; |
---|
| 781 | vtable(ivec)=vector(2)+shifty; |
---|
| 782 | end |
---|
| 783 | result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output |
---|
| 784 | |
---|
| 785 | |
---|
| 786 | function [vector,F] = SUBPIXGAUSS (result_conv,x,y) |
---|
| 787 | vector=[0 0]; %default |
---|
| 788 | F=0; |
---|
| 789 | [npy,npx]=size(result_conv); |
---|
| 790 | |
---|
| 791 | % if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1) |
---|
| 792 | %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion Israel Institute of Technology |
---|
| 793 | %http://urapiv.wordpress.com |
---|
| 794 | peaky = y; |
---|
| 795 | if y <= npy-1 && y >= 1 |
---|
| 796 | f0 = log(result_conv(y,x)); |
---|
| 797 | f1 = real(log(result_conv(y-1,x))); |
---|
| 798 | f2 = real(log(result_conv(y+1,x))); |
---|
| 799 | peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2); |
---|
| 800 | else |
---|
| 801 | F=-2; % warning flag for vector truncated by the limited search box |
---|
| 802 | end |
---|
| 803 | peakx=x; |
---|
| 804 | if x <= npx-1 && x >= 1 |
---|
| 805 | f0 = log(result_conv(y,x)); |
---|
| 806 | f1 = real(log(result_conv(y,x-1))); |
---|
| 807 | f2 = real(log(result_conv(y,x+1))); |
---|
| 808 | peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2); |
---|
| 809 | else |
---|
| 810 | F=-2; % warning flag for vector truncated by the limited search box |
---|
| 811 | end |
---|
| 812 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|
| 813 | % else |
---|
| 814 | % vector=[NaN NaN]; |
---|
| 815 | % end |
---|
| 816 | |
---|
| 817 | function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y) |
---|
| 818 | vector=[0 0]; %default |
---|
| 819 | F=-2; |
---|
| 820 | peaky=y; |
---|
| 821 | peakx=x; |
---|
| 822 | [npy,npx]=size(result_conv); |
---|
| 823 | if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1) |
---|
| 824 | F=0; |
---|
| 825 | for i=-1:1 |
---|
| 826 | for j=-1:1 |
---|
| 827 | %following 15 lines based on |
---|
| 828 | %H. Nobach Æ M. Honkanen (2005) |
---|
| 829 | %Two-dimensional Gaussian regression for sub-pixel displacement |
---|
| 830 | %estimation in particle image velocimetry or particle position |
---|
| 831 | %estimation in particle tracking velocimetry |
---|
| 832 | %Experiments in Fluids (2005) 38: 511515 |
---|
| 833 | c10(j+2,i+2)=i*log(result_conv(y+j, x+i)); |
---|
| 834 | c01(j+2,i+2)=j*log(result_conv(y+j, x+i)); |
---|
| 835 | c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i)); |
---|
| 836 | c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i)); |
---|
| 837 | c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i)); |
---|
| 838 | %c00(j+2,i+2)=(5-3*i^2-3*j^2)*log(result_conv_norm(maxY+j, maxX+i)); |
---|
| 839 | end |
---|
| 840 | end |
---|
| 841 | c10=(1/6)*sum(sum(c10)); |
---|
| 842 | c01=(1/6)*sum(sum(c01)); |
---|
| 843 | c11=(1/4)*sum(sum(c11)); |
---|
| 844 | c20=(1/6)*sum(sum(c20)); |
---|
| 845 | c02=(1/6)*sum(sum(c02)); |
---|
| 846 | deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2); |
---|
| 847 | deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2); |
---|
| 848 | if abs(deltax)<1 |
---|
| 849 | peakx=x+deltax; |
---|
| 850 | end |
---|
| 851 | if abs(deltay)<1 |
---|
| 852 | peaky=y+deltay; |
---|
| 853 | end |
---|
| 854 | end |
---|
| 855 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1]; |
---|
| 856 | |
---|
| 857 | |
---|