[242] | 1 | % To develop.... |
---|
| 2 | function [Data,errormsg]= civ_uvmat(Param,ncfile) |
---|
[246] | 3 | errormsg=''; |
---|
[242] | 4 | Data.ListGlobalAttribute={'Conventions','Program','CivStage'}; |
---|
| 5 | Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes |
---|
| 6 | Data.Program='civ_uvmat'; |
---|
| 7 | Data.CivStage=0;%default |
---|
[246] | 8 | ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'}; |
---|
| 9 | ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'}; |
---|
[248] | 10 | mask=''; |
---|
| 11 | maskname='';%default |
---|
[251] | 12 | test_civx=0;%default |
---|
[242] | 13 | |
---|
| 14 | %% Civ1 |
---|
| 15 | if isfield (Param,'Civ1') |
---|
| 16 | par_civ1=Param.Civ1; |
---|
| 17 | image1=imread(par_civ1.filename_ima_a); |
---|
| 18 | image2=imread(par_civ1.filename_ima_b); |
---|
| 19 | stepx=str2num(par_civ1.dx); |
---|
| 20 | stepy=str2num(par_civ1.dy); |
---|
| 21 | ibx2=ceil(str2num(par_civ1.ibx)/2); |
---|
| 22 | iby2=ceil(str2num(par_civ1.iby)/2); |
---|
| 23 | isx2=ceil(str2num(par_civ1.isx)/2); |
---|
| 24 | isy2=ceil(str2num(par_civ1.isy)/2); |
---|
| 25 | shiftx=str2num(par_civ1.shiftx); |
---|
| 26 | shifty=str2num(par_civ1.shifty); |
---|
[246] | 27 | miniy=max(1+isy2+shifty,1+iby2); |
---|
[242] | 28 | minix=max(1+isx2-shiftx,1+ibx2); |
---|
[246] | 29 | maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2); |
---|
[242] | 30 | maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2); |
---|
| 31 | [GridX,GridY]=meshgrid(minix:stepx:maxix,miniy:stepy:maxiy); |
---|
| 32 | PointCoord(:,1)=reshape(GridX,[],1); |
---|
| 33 | PointCoord(:,2)=reshape(GridY,[],1); |
---|
[252] | 34 | if isfield(par_civ1,'maskname') && ~isempty(par_civ1.maskname) |
---|
[248] | 35 | maskname=par_civ1.maskname; |
---|
| 36 | mask=imread(maskname); |
---|
| 37 | end |
---|
[242] | 38 | % caluclate velocity data (y and v in indices, reverse to y component) |
---|
[248] | 39 | [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho), mask); |
---|
[242] | 40 | list_param=(fieldnames(par_civ1))'; |
---|
| 41 | list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'}; |
---|
| 42 | index=zeros(size(list_param)); |
---|
| 43 | for ilist=1:length(list_remove) |
---|
| 44 | index=strcmp(list_remove{ilist},list_param); |
---|
| 45 | if ~isempty(find(index,1)) |
---|
| 46 | list_param(index)=[]; |
---|
| 47 | end |
---|
| 48 | end |
---|
| 49 | for ilist=1:length(list_param) |
---|
| 50 | Civ1_param{ilist}=['Civ1_' list_param{ilist}]; |
---|
| 51 | eval(['Data.Civ1_' list_param{ilist} '=Param.Civ1.' list_param{ilist} ';']) |
---|
| 52 | end |
---|
| 53 | if isfield(Data,'Civ1_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile') |
---|
| 54 | Data.Civ1_gridname=''; |
---|
| 55 | end |
---|
| 56 | if isfield(Data,'Civ1_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile') |
---|
| 57 | Data.Civ1_maskname=''; |
---|
| 58 | end |
---|
| 59 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param {'Civ1_Time','Civ1_Dt'}]; |
---|
| 60 | Data.Civ1_Time=str2double(par_civ1.T0); |
---|
| 61 | Data.Civ1_Dt=str2double(par_civ1.Dt); |
---|
| 62 | Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};% cell array containing the names of the fields to record |
---|
| 63 | Data.VarDimName={'nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'}; |
---|
| 64 | Data.VarAttribute{1}.Role='coord_x'; |
---|
| 65 | Data.VarAttribute{2}.Role='coord_y'; |
---|
| 66 | Data.VarAttribute{3}.Role='vector_x'; |
---|
| 67 | Data.VarAttribute{4}.Role='vector_y'; |
---|
| 68 | Data.VarAttribute{5}.Role='warnflag'; |
---|
| 69 | Data.Civ1_X=reshape(xtable,[],1); |
---|
| 70 | Data.Civ1_Y=reshape(size(image1,1)-ytable+1,[],1); |
---|
| 71 | Data.Civ1_U=reshape(utable,[],1); |
---|
| 72 | Data.Civ1_V=reshape(-vtable,[],1); |
---|
| 73 | Data.Civ1_C=reshape(ctable,[],1); |
---|
| 74 | Data.Civ1_F=reshape(F,[],1); |
---|
| 75 | Data.CivStage=1; |
---|
| 76 | else |
---|
[251] | 77 | Data=nc2struct(ncfile,'ListGlobalAttribute','absolut_time_T0'); |
---|
| 78 | |
---|
| 79 | % read Civx data |
---|
| 80 | if ~isempty(Data.absolut_time_T0')%read civx file |
---|
| 81 | % var={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF';... |
---|
| 82 | % 'vec_X','vec_Y','vec_U','vec_V','vec_C','vec_F','vec_FixFlag'}; |
---|
| 83 | |
---|
| 84 | %var=varcivx_generator('velocity','Civ1');%determine the names of constants and variables to read |
---|
| 85 | test_civx=1; |
---|
| 86 | [Data,vardetect,ichoice]=nc2struct(ncfile);%read the variables in the netcdf file |
---|
| 87 | % Data.ListGlobalAttribute=[{'Conventions','Program','CivStage'} Data.ListGlobalAttribute {'Civ1_Time','Civ1_Dt'}]; |
---|
| 88 | % Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes |
---|
| 89 | % Data.Program='civ_uvmat'; |
---|
| 90 | % Data.Civ1_Time=double(Data.absolut_time_T0); |
---|
| 91 | % Data.Civ1_Dt=double(Data.dt); |
---|
| 92 | % Data.VarDimName={'nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'}; |
---|
| 93 | % Data.VarAttribute{1}.Role='coord_x'; |
---|
| 94 | % Data.VarAttribute{2}.Role='coord_y'; |
---|
| 95 | % Data.VarAttribute{3}.Role='vector_x'; |
---|
| 96 | % Data.VarAttribute{4}.Role='vector_y'; |
---|
| 97 | % Data.VarAttribute{5}.Role='ancillary'; |
---|
| 98 | % Data.VarAttribute{6}.Role='warnflag'; |
---|
| 99 | % Data.VarAttribute{7}.Role='errorflag'; |
---|
| 100 | % Data.CivStage=1; |
---|
[246] | 101 | else |
---|
[251] | 102 | |
---|
| 103 | if isfield(Param,'Fix1') |
---|
| 104 | Data=nc2struct(ncfile,ListVarCiv1);%read civ1 data in the existing netcdf file |
---|
| 105 | else |
---|
| 106 | Data=nc2struct(ncfile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file |
---|
| 107 | end |
---|
[246] | 108 | end |
---|
| 109 | if isfield(Data,'Txt') |
---|
| 110 | msgbox_uvmat('ERROR',Data.Txt) |
---|
| 111 | return |
---|
| 112 | end |
---|
[242] | 113 | end |
---|
| 114 | |
---|
| 115 | %% Fix1 |
---|
| 116 | if isfield (Param,'Fix1') |
---|
| 117 | ListFixParam=fieldnames(Param.Fix1); |
---|
| 118 | for ilist=1:length(ListFixParam) |
---|
| 119 | ParamName=ListFixParam{ilist}; |
---|
| 120 | ListName=['Fix1_' ParamName]; |
---|
| 121 | eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']) |
---|
| 122 | eval(['Data.' ListName '=Param.Fix1.' ParamName ';']) |
---|
| 123 | end |
---|
| 124 | % Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Fix1_WarnFlags','Fix1_TreshCorr','Fix1_TreshVel','Fix1_UpperBoundTest'}]; |
---|
| 125 | % Data.Fix1_WarnFlags=Param.Fix1.WarnFlags; |
---|
| 126 | % Data.Fix1_ThreshCorr=Param.Fix1.ThreshCorr; |
---|
| 127 | % Data.Fix1_ThreshVel=Param.Fix1.ThreshVel; |
---|
| 128 | % Data.Fix1_UpperBoundTest=Param.Fix1.UpperBoundTest; |
---|
[251] | 129 | |
---|
| 130 | if test_civx |
---|
| 131 | if ~isfield(Data,'fix') |
---|
| 132 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix']; |
---|
| 133 | Data.fix=1; |
---|
| 134 | Data.ListVarName=[Data.ListVarName {'vec_FixFlag'}]; |
---|
[252] | 135 | Data.VarDimName=[Data.VarDimName {'nb_vectors'}]; |
---|
[251] | 136 | end |
---|
| 137 | Data.vec_FixFlag=fix_uvmat(Param.Fix1,Data.vec_F,Data.vec_C,Data.vec_U,Data.vec_V,Data.vec_X,Data.vec_Y); |
---|
| 138 | else |
---|
| 139 | Data.ListVarName=[Data.ListVarName {'Civ1_FF'}]; |
---|
| 140 | Data.VarDimName=[Data.VarDimName {'nbvec1'}]; |
---|
| 141 | nbvar=length(Data.ListVarName); |
---|
| 142 | Data.VarAttribute{nbvar}.Role='errorflag'; |
---|
| 143 | Data.Civ1_FF=fix_uvmat(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V); |
---|
| 144 | Data.CivStage=2; |
---|
| 145 | end |
---|
[242] | 146 | end |
---|
| 147 | %% Patch1 |
---|
| 148 | if isfield (Param,'Patch1') |
---|
| 149 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}]; |
---|
| 150 | Data.Patch1_Rho=str2double(Param.Patch1.Rho); |
---|
| 151 | Data.Patch1_Threshold=str2double(Param.Patch1.Threshold); |
---|
| 152 | Data.Patch1_SubDomain=str2double(Param.Patch1.SubDomain); |
---|
[246] | 153 | 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'}]; |
---|
| 154 | Data.VarDimName=[Data.VarDimName {'NbVec1','NbVec1',{'NbSubDomain1','Two'},{'NbSubDomain1','Two'},'NbSubDomain1',... |
---|
| 155 | {'NbVec1Sub','NbSubDomain1'},{'NbVec1Sub','NbSubDomain1'},{'Nbtps1','NbSubDomain1'},{'Nbtps1','NbSubDomain1'}}]; |
---|
[242] | 156 | nbvar=length(Data.ListVarName); |
---|
| 157 | Data.VarAttribute{nbvar-1}.Role='vector_x'; |
---|
| 158 | Data.VarAttribute{nbvar}.Role='vector_y'; |
---|
[246] | 159 | Data.Civ1_U_Diff=zeros(size(Data.Civ1_X)); |
---|
| 160 | Data.Civ1_V_Diff=zeros(size(Data.Civ1_X)); |
---|
[242] | 161 | if isfield(Data,'Civ1_FF') |
---|
| 162 | ind_good=find(Data.Civ1_FF==0); |
---|
| 163 | else |
---|
| 164 | ind_good=1:numel(Data.Civ1_X); |
---|
| 165 | end |
---|
[246] | 166 | [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]=... |
---|
[242] | 167 | patch_uvmat(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] | 168 | Data.Civ1_U_Diff(ind_good)=Data.Civ1_U(ind_good)-Ures; |
---|
| 169 | Data.Civ1_V_Diff(ind_good)=Data.Civ1_V(ind_good)-Vres; |
---|
| 170 | Data.Civ1_FF(ind_good)=FFres; |
---|
| 171 | Data.CivStage=3; |
---|
[242] | 172 | end |
---|
[248] | 173 | |
---|
| 174 | %% Civ2 |
---|
| 175 | if isfield (Param,'Civ2') |
---|
[252] | 176 | par_civ2=Param.Civ2; |
---|
[248] | 177 | image1=imread(par_civ2.filename_ima_a); |
---|
| 178 | image2=imread(par_civ2.filename_ima_b); |
---|
| 179 | stepx=str2num(par_civ2.dx); |
---|
| 180 | stepy=str2num(par_civ2.dy); |
---|
| 181 | ibx2=ceil(str2num(par_civ2.ibx)/2); |
---|
| 182 | iby2=ceil(str2num(par_civ2.iby)/2); |
---|
| 183 | isx2=4; |
---|
| 184 | isy2=4; |
---|
| 185 | % shiftx=str2num(par_civ1.shiftx); |
---|
| 186 | % shifty=str2num(par_civ1.shifty); |
---|
| 187 | % TO GET shift from par_civ2.filename_nc1 |
---|
| 188 | miniy=max(1+isy2+shifty,1+iby2); |
---|
| 189 | minix=max(1+isx2-shiftx,1+ibx2); |
---|
| 190 | maxiy=min(size(image1,1)-isy2+shifty,size(image1,1)-iby2); |
---|
| 191 | maxix=min(size(image1,2)-isx2-shiftx,size(image1,2)-ibx2); |
---|
| 192 | [GridX,GridY]=meshgrid(minix:stepx:maxix,miniy:stepy:maxiy); |
---|
| 193 | PointCoord(:,1)=reshape(GridX,[],1); |
---|
| 194 | PointCoord(:,2)=reshape(GridY,[],1); |
---|
| 195 | if ~isempty(par_civ2.maskname)&& ~strcmp(maskname,par_civ2.maskname)% mask exist, not already read in civ1 |
---|
| 196 | mask=imread(par_civ2.maskname); |
---|
| 197 | end |
---|
| 198 | % caluclate velocity data (y and v in indices, reverse to y component) |
---|
| 199 | [xtable ytable utable vtable ctable F] = pivlab (image1,image2,ibx2,iby2,isx2,isy2,shiftx,-shifty,PointCoord,str2num(par_civ1.rho),mask); |
---|
| 200 | list_param=(fieldnames(par_civ1))'; |
---|
| 201 | list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'}; |
---|
| 202 | index=zeros(size(list_param)); |
---|
| 203 | for ilist=1:length(list_remove) |
---|
| 204 | index=strcmp(list_remove{ilist},list_param); |
---|
| 205 | if ~isempty(find(index,1)) |
---|
| 206 | list_param(index)=[]; |
---|
| 207 | end |
---|
| 208 | end |
---|
| 209 | for ilist=1:length(list_param) |
---|
| 210 | Civ1_param{ilist}=['Civ1_' list_param{ilist}]; |
---|
| 211 | eval(['Data.Civ1_' list_param{ilist} '=Param.Civ1.' list_param{ilist} ';']) |
---|
| 212 | end |
---|
| 213 | if isfield(Data,'Civ1_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile') |
---|
| 214 | Data.Civ1_gridname=''; |
---|
| 215 | end |
---|
| 216 | if isfield(Data,'Civ1_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile') |
---|
| 217 | Data.Civ1_maskname=''; |
---|
| 218 | end |
---|
| 219 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param {'Civ1_Time','Civ1_Dt'}]; |
---|
| 220 | Data.Civ1_Time=str2double(par_civ1.T0); |
---|
| 221 | Data.Civ1_Dt=str2double(par_civ1.Dt); |
---|
| 222 | Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};% cell array containing the names of the fields to record |
---|
| 223 | Data.VarDimName={'nbvec1','nbvec1','nbvec1','nbvec1','nbvec1','nbvec1'}; |
---|
| 224 | Data.VarAttribute{1}.Role='coord_x'; |
---|
| 225 | Data.VarAttribute{2}.Role='coord_y'; |
---|
| 226 | Data.VarAttribute{3}.Role='vector_x'; |
---|
| 227 | Data.VarAttribute{4}.Role='vector_y'; |
---|
| 228 | Data.VarAttribute{5}.Role='warnflag'; |
---|
| 229 | Data.Civ1_X=reshape(xtable,[],1); |
---|
| 230 | Data.Civ1_Y=reshape(size(image1,1)-ytable+1,[],1); |
---|
| 231 | Data.Civ1_U=reshape(utable,[],1); |
---|
| 232 | Data.Civ1_V=reshape(-vtable,[],1); |
---|
| 233 | Data.Civ1_C=reshape(ctable,[],1); |
---|
| 234 | Data.Civ1_F=reshape(F,[],1); |
---|
| 235 | Data.CivStage=Data.CivStage+1; |
---|
| 236 | end |
---|
[252] | 237 | |
---|
| 238 | %% Fix2 |
---|
| 239 | if isfield (Param,'Fix2') |
---|
| 240 | ListFixParam=fieldnames(Param.Fix2); |
---|
| 241 | for ilist=1:length(ListFixParam) |
---|
| 242 | ParamName=ListFixParam{ilist}; |
---|
| 243 | ListName=['Fix1_' ParamName]; |
---|
| 244 | eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];']) |
---|
| 245 | eval(['Data.' ListName '=Param.Fix2.' ParamName ';']) |
---|
| 246 | end |
---|
| 247 | if test_civx |
---|
| 248 | if ~isfield(Data,'fix2') |
---|
| 249 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2']; |
---|
| 250 | Data.fix2=1; |
---|
| 251 | Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}]; |
---|
| 252 | Data.VarDimName=[Data.VarDimName {'nb_vectors2'}]; |
---|
| 253 | end |
---|
| 254 | Data.vec_FixFlag=fix_uvmat(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y); |
---|
| 255 | else |
---|
| 256 | Data.ListVarName=[Data.ListVarName {'Civ2_FF'}]; |
---|
| 257 | Data.VarDimName=[Data.VarDimName {'nbvec2'}]; |
---|
| 258 | nbvar=length(Data.ListVarName); |
---|
| 259 | Data.VarAttribute{nbvar}.Role='errorflag'; |
---|
| 260 | Data.Civ2_FF=fix_uvmat(Param.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V); |
---|
| 261 | Data.CivStage=5; |
---|
| 262 | end |
---|
| 263 | end |
---|
| 264 | |
---|
| 265 | %% Patch2 |
---|
| 266 | if isfield (Param,'Patch2') |
---|
| 267 | Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch2_Rho','Patch2_Threshold','Patch2_SubDomain'}]; |
---|
| 268 | Data.Patch2_Rho=str2double(Param.Patch2.Rho); |
---|
| 269 | Data.Patch2_Threshold=str2double(Param.Patch2.Threshold); |
---|
| 270 | Data.Patch2_SubDomain=str2double(Param.Patch2.SubDomain); |
---|
| 271 | 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'}]; |
---|
| 272 | Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2',{'NbSubDomain2','Two'},{'NbSubDomain2','Two'},'NbSubDomain2',... |
---|
| 273 | {'NbVec2Sub','NbSubDomain2'},{'NbVec2Sub','NbSubDomain2'},{'Nbtps2','NbSubDomain2'},{'Nbtps2','NbSubDomain2'}}]; |
---|
| 274 | nbvar=length(Data.ListVarName); |
---|
| 275 | Data.VarAttribute{nbvar-1}.Role='vector_x'; |
---|
| 276 | Data.VarAttribute{nbvar}.Role='vector_y'; |
---|
| 277 | Data.Civ2_U_Diff=zeros(size(Data.Civ2_X)); |
---|
| 278 | Data.Civ2_V_Diff=zeros(size(Data.Civ2_X)); |
---|
| 279 | if isfield(Data,'Civ2_FF') |
---|
| 280 | ind_good=find(Data.Civ2_FF==0); |
---|
| 281 | else |
---|
| 282 | ind_good=1:numel(Data.Civ2_X); |
---|
| 283 | end |
---|
| 284 | [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]=... |
---|
| 285 | patch_uvmat(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); |
---|
| 286 | Data.Civ2_U_Diff(ind_good)=Data.Civ2_U(ind_good)-Ures; |
---|
| 287 | Data.Civ2_V_Diff(ind_good)=Data.Civ2_V(ind_good)-Vres; |
---|
| 288 | Data.Civ2_FF(ind_good)=FFres; |
---|
| 289 | Data.CivStage=3; |
---|
| 290 | end |
---|
| 291 | |
---|
[242] | 292 | %% write result |
---|
[246] | 293 | % 'TESTcalc' |
---|
| 294 | % [DataOut,errormsg]=calc_field('velocity',Data) |
---|
| 295 | if exist('ncfile','var') |
---|
[242] | 296 | errormsg=struct2nc(ncfile,Data); |
---|
[246] | 297 | end |
---|
[242] | 298 | |
---|
| 299 | |
---|
| 300 | %'RUN_FIX': function for fixing velocity fields: |
---|
| 301 | %----------------------------------------------- |
---|
| 302 | % RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref) |
---|
| 303 | % |
---|
| 304 | %filename: name of the netcdf file (used as input and output) |
---|
| 305 | %field: structure specifying the names of the fields to fix (depending on civ1 or civ2) |
---|
| 306 | %.vel_type='civ1' or 'civ2'; |
---|
| 307 | %.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1); |
---|
| 308 | %.fixflag=name of fix flag variable ('vec_FixFlag' for civ1) |
---|
| 309 | %flagindex: flag specifying which values of vec_f are removed: |
---|
| 310 | % if flagindex(1)=1: vec_f=-2 vectors are removed |
---|
| 311 | % if flagindex(2)=1: vec_f=3 vectors are removed |
---|
| 312 | % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2) |
---|
| 313 | %iter=1 for civ1 fields and iter=2 for civ2 fields |
---|
| 314 | %thresh_vecC: threshold in the image correlation vec_C |
---|
| 315 | %flag_mask: =1 mask used to remove vectors (0 else) |
---|
| 316 | %maskname: name of the mask image file for fix |
---|
| 317 | %thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists |
---|
| 318 | %inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold |
---|
| 319 | %fileref: .nc file name for a reference velocity (='': refrence 0 used) |
---|
| 320 | %fieldref: 'civ1','filter1'...feld used in fileref |
---|
| 321 | |
---|
[251] | 322 | function FF=fix_uvmat(Param,F,C,U,V,X,Y) |
---|
[242] | 323 | FF=zeros(size(F));%default |
---|
| 324 | |
---|
| 325 | %criterium on warn flags |
---|
| 326 | if isfield (Param,'WarnFlags') |
---|
| 327 | for iflag=1:numel(Param.WarnFlags) |
---|
| 328 | FF=(FF==1| F==Param.WarnFlags(iflag)); |
---|
| 329 | end |
---|
| 330 | end |
---|
| 331 | |
---|
| 332 | %criterium on correlation values |
---|
| 333 | if isfield (Param,'LowerBoundCorr') |
---|
| 334 | FF=FF==1 | C<Param.LowerBoundCorr; |
---|
| 335 | end |
---|
| 336 | |
---|
| 337 | if isfield (Param,'LowerBoundVel')&& ~isequal(Param.LowerBoundVel,0) |
---|
| 338 | thresh=Param.LowerBoundVel*Param.LowerBoundVel; |
---|
| 339 | FF=FF==1 | (U.*U+V.*V)<thresh; |
---|
| 340 | end |
---|
| 341 | if isfield (Param,'UpperBoundVel')&& ~isequal(Param.UpperBoundVel,0) |
---|
| 342 | thresh=Param.UpperBoundVel*Param.UpperBoundVel; |
---|
| 343 | FF=FF==1 | (U.*U+V.*V)>thresh; |
---|
| 344 | end |
---|
[251] | 345 | if isfield(Param,'MaskName') |
---|
| 346 | M=imread(Param.MaskName); |
---|
| 347 | nxy=size(M); |
---|
| 348 | M=reshape(M,1,[]); |
---|
| 349 | rangx0=[0.5 nxy(2)-0.5]; |
---|
| 350 | rangy0=[0.5 nxy(1)-0.5]; |
---|
| 351 | vec_x1=X-U/2;%beginning points |
---|
| 352 | vec_x2=X+U/2;%end points of vectors |
---|
| 353 | vec_y1=Y-V/2;%beginning points |
---|
| 354 | vec_y2=Y+V/2;%end points of vectors |
---|
| 355 | indx=1+round((nxy(2)-1)*(vec_x1-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x1 |
---|
| 356 | indy=1+round((nxy(1)-1)*(vec_y1-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y1 |
---|
| 357 | test_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside |
---|
| 358 | indx=indx.*test_in+(1-test_in); %replace indx by 1 out of the mask range |
---|
| 359 | indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range |
---|
| 360 | ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector |
---|
| 361 | Mvalues=M(ICOMB); |
---|
| 362 | flag7b=((20 < Mvalues) & (Mvalues < 200))| ~test_in'; |
---|
| 363 | indx=1+round((nxy(2)-1)*(vec_x2-rangx0(1))/(rangx0(2)-rangx0(1)));% image index x at abcissa vec_x2 |
---|
| 364 | indy=1+round((nxy(1)-1)*(vec_y2-rangy0(1))/(rangy0(2)-rangy0(1)));% image index y at ordinate vec_y2 |
---|
| 365 | test_in=~(indx < 1 |indy < 1 | indx > nxy(2) |indy > nxy(1)); %=0 out of the mask image, 1 inside |
---|
| 366 | indx=indx.*test_in+(1-test_in); %replace indx by 1 out of the mask range |
---|
| 367 | indy=indy.*test_in+(1-test_in); %replace indy by 1 out of the mask range |
---|
| 368 | ICOMB=((indx-1)*nxy(1)+(nxy(1)+1-indy));%determine the indices in the image reshaped in a Matlab vector |
---|
| 369 | Mvalues=M(ICOMB); |
---|
| 370 | flag7e=((Mvalues > 20) & (Mvalues < 200))| ~test_in'; |
---|
| 371 | FF=FF==1 |(flag7b|flag7e)'; |
---|
| 372 | end |
---|
| 373 | % flag7=0; |
---|
| 374 | % end |
---|
| 375 | |
---|
| 376 | |
---|
[246] | 377 | FF=double(FF); |
---|
[242] | 378 | % |
---|
| 379 | % % criterium on velocity values |
---|
| 380 | % delta_u=Field.U;%default without ref file |
---|
| 381 | % delta_v=Field.V; |
---|
| 382 | % if exist('fileref','var') && ~isempty(fileref) |
---|
| 383 | % if ~exist(fileref,'file') |
---|
| 384 | % error='reference file not found in RUN_FIX.m'; |
---|
| 385 | % display(error); |
---|
| 386 | % return |
---|
| 387 | % end |
---|
| 388 | % FieldRef=read_civxdata(fileref,[],fieldref); |
---|
| 389 | % if isfield(FieldRef,'FF') |
---|
| 390 | % index_true=find(FieldRef.FF==0); |
---|
| 391 | % FieldRef.X=FieldRef.X(index_true); |
---|
| 392 | % FieldRef.Y=FieldRef.Y(index_true); |
---|
| 393 | % FieldRef.U=FieldRef.U(index_true); |
---|
| 394 | % FieldRef.V=FieldRef.V(index_true); |
---|
| 395 | % end |
---|
| 396 | % if ~isfield(FieldRef,'X') || ~isfield(FieldRef,'Y') || ~isfield(FieldRef,'U') || ~isfield(FieldRef,'V') |
---|
| 397 | % error='reference file is not a velocity field in RUN_FIX.m '; %bad input file |
---|
| 398 | % return |
---|
| 399 | % end |
---|
| 400 | % if length(FieldRef.X)<=1 |
---|
| 401 | % errordlg('reference field with one vector or less in RUN_FIX.m') |
---|
| 402 | % return |
---|
| 403 | % end |
---|
| 404 | % vec_U_ref=griddata_uvmat(FieldRef.X,FieldRef.Y,FieldRef.U,Field.X,Field.Y); %interpolate vectors in the ref field |
---|
| 405 | % 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 |
---|
| 406 | % delta_u=Field.U-vec_U_ref;%take the difference with the interpolated ref field |
---|
| 407 | % delta_v=Field.V-vec_V_ref; |
---|
| 408 | % end |
---|
| 409 | % thresh_vel_x=thresh_vel; |
---|
| 410 | % thresh_vel_y=thresh_vel; |
---|
| 411 | % if isequal(inf_sup,1) |
---|
| 412 | % flag5=abs(delta_u)<thresh_vel_x & abs(delta_v)<thresh_vel_y &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1); |
---|
| 413 | % elseif isequal(inf_sup,2) |
---|
| 414 | % flag5=(abs(delta_u)>thresh_vel_x | abs(delta_v)>thresh_vel_y) &(flag1~=1)&(flag2~=1)&(flag3~=1)&(flag4~=1); |
---|
| 415 | % end |
---|
| 416 | % |
---|
| 417 | % % flag7 introduce a grey mask, matrix M |
---|
[251] | 418 | |
---|
[242] | 419 | % flagmagenta=flag1|flag2|flag3|flag4|flag5|flag7; |
---|
| 420 | % fixflag_unit=Field.FF-10*floor(Field.FF/10); %unity term of fix_flag |
---|
| 421 | |
---|
| 422 | |
---|
| 423 | |
---|
| 424 | %------------------------------------------------------------------------ |
---|
| 425 | % patch function |
---|
[246] | 426 | function [SubRangx,SubRangy,nbpoints,FF,U_smooth,V_smooth,X_tps,Y_tps,U_tps,V_tps] =patch_uvmat(X,Y,U,V,Rho,Threshold,SubDomain) |
---|
[242] | 427 | %subdomain decomposition |
---|
| 428 | warning off |
---|
[243] | 429 | U=reshape(U,[],1); |
---|
| 430 | V=reshape(V,[],1); |
---|
| 431 | X=reshape(X,[],1); |
---|
| 432 | Y=reshape(Y,[],1); |
---|
[242] | 433 | nbvec=numel(X); |
---|
[246] | 434 | NbSubDomain=ceil(nbvec/SubDomain); |
---|
| 435 | MinX=min(X); |
---|
| 436 | MinY=min(Y); |
---|
| 437 | MaxX=max(X); |
---|
| 438 | MaxY=max(Y); |
---|
[242] | 439 | RangX=MaxX-MinX; |
---|
| 440 | RangY=MaxY-MinY; |
---|
[246] | 441 | AspectRatio=RangY/RangX; |
---|
| 442 | NbSubDomainX=max(floor(sqrt(NbSubDomain/AspectRatio)),1); |
---|
| 443 | NbSubDomainY=max(floor(sqrt(NbSubDomain*AspectRatio)),1); |
---|
| 444 | NbSubDomain=NbSubDomainX*NbSubDomainY; |
---|
[242] | 445 | SizX=RangX/NbSubDomainX;%width of subdomains |
---|
| 446 | SizY=RangY/NbSubDomainY;%height of subdomains |
---|
| 447 | CentreX=linspace(MinX+SizX/2,MaxX-SizX/2,NbSubDomainX); |
---|
| 448 | CentreY=linspace(MinY+SizY/2,MaxY-SizY/2,NbSubDomainY); |
---|
[246] | 449 | [CentreX,CentreY]=meshgrid(CentreX,CentreY); |
---|
| 450 | CentreY=reshape(CentreY,1,[]); |
---|
| 451 | CentreX=reshape(CentreX,1,[]); |
---|
| 452 | 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) |
---|
| 453 | U_tps_sub=zeros(length(X),NbSubDomain);%default spline |
---|
| 454 | V_tps_sub=zeros(length(X),NbSubDomain);%default spline |
---|
| 455 | U_smooth=zeros(length(X),1); |
---|
| 456 | V_smooth=zeros(length(X),1); |
---|
| 457 | |
---|
| 458 | nb_select=zeros(length(X),1); |
---|
[243] | 459 | FF=zeros(length(X),1); |
---|
[246] | 460 | test_empty=zeros(1,NbSubDomain); |
---|
| 461 | for isub=1:NbSubDomain |
---|
| 462 | SubRangx(isub,:)=[CentreX(isub)-SizX/2 CentreX(isub)+SizX/2]; |
---|
| 463 | SubRangy(isub,:)=[CentreY(isub)-SizY/2 CentreY(isub)+SizY/2]; |
---|
| 464 | ind_sel_previous=[]; |
---|
| 465 | ind_sel=0; |
---|
| 466 | while numel(ind_sel)>numel(ind_sel_previous) %increase the subdomain during four iterations at most |
---|
| 467 | ind_sel_previous=ind_sel; |
---|
| 468 | ind_sel=find(X>SubRangx(isub,1) & X<SubRangx(isub,2) & Y>SubRangy(isub,1) & Y<SubRangy(isub,2)); |
---|
| 469 | % if no vector in the subdomain, skip the subdomain |
---|
| 470 | if isempty(ind_sel) |
---|
| 471 | test_empty(isub)=1; |
---|
| 472 | U_tps(1,isub)=0;%define U_tps and V_tps by default |
---|
| 473 | V_tps(1,isub)=0; |
---|
| 474 | break |
---|
| 475 | % if too few selected vectors, increase the subrange for next iteration |
---|
| 476 | elseif numel(ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); |
---|
| 477 | SubRangx(isub,1)=SubRangx(isub,1)-SizX/4; |
---|
| 478 | SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; |
---|
| 479 | SubRangy(isub,1)=SubRangy(isub,1)-SizY/4; |
---|
| 480 | SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; |
---|
| 481 | else |
---|
| 482 | [U_smooth_sub,U_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),U(ind_sel),rho); |
---|
| 483 | [V_smooth_sub,V_tps_sub]=tps_coeff(X(ind_sel),Y(ind_sel),V(ind_sel),rho); |
---|
| 484 | UDiff=U_smooth_sub-U(ind_sel); |
---|
| 485 | VDiff=V_smooth_sub-V(ind_sel); |
---|
| 486 | NormDiff=UDiff.*UDiff+VDiff.*VDiff; |
---|
| 487 | FF(ind_sel)=20*(NormDiff>Threshold);%put FF value to 20 to identify the criterium of elimmination |
---|
| 488 | ind_ind_sel=find(FF(ind_sel)==0); % select the indices of ind_sel corresponding to the remaining vectors |
---|
| 489 | % no value exceeds threshold, the result is recorded |
---|
| 490 | if isequal(numel(ind_ind_sel),numel(ind_sel)) |
---|
| 491 | U_smooth(ind_sel)=U_smooth(ind_sel)+U_smooth_sub; |
---|
| 492 | V_smooth(ind_sel)=V_smooth(ind_sel)+V_smooth_sub; |
---|
| 493 | nbpoints(isub)=numel(ind_sel); |
---|
| 494 | X_tps(1:nbpoints(isub),isub)=X(ind_sel); |
---|
| 495 | Y_tps(1:nbpoints(isub),isub)=Y(ind_sel); |
---|
| 496 | U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; |
---|
| 497 | V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; |
---|
| 498 | nb_select(ind_sel)=nb_select(ind_sel)+1; |
---|
| 499 | display('good') |
---|
| 500 | break |
---|
| 501 | % too few selected vectors, increase the subrange for next iteration |
---|
| 502 | elseif numel(ind_ind_sel)<SubDomain/4 && ~isequal( ind_sel,ind_sel_previous); |
---|
| 503 | SubRangx(isub,1)=SubRangx(isub,1)-SizX/4; |
---|
| 504 | SubRangx(isub,2)=SubRangx(isub,2)+SizX/4; |
---|
| 505 | SubRangy(isub,1)=SubRangy(isub,1)-SizY/4; |
---|
| 506 | SubRangy(isub,2)=SubRangy(isub,2)+SizY/4; |
---|
| 507 | % display('fewsmooth') |
---|
| 508 | % interpolation-smoothing is done again with the selected vectors |
---|
[242] | 509 | else |
---|
[246] | 510 | [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); |
---|
| 511 | [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); |
---|
| 512 | U_smooth(ind_sel(ind_ind_sel))=U_smooth(ind_sel(ind_ind_sel))+U_smooth_sub; |
---|
| 513 | V_smooth(ind_sel(ind_ind_sel))=V_smooth(ind_sel(ind_ind_sel))+V_smooth_sub; |
---|
| 514 | nbpoints(isub)=numel(ind_ind_sel); |
---|
| 515 | X_tps(1:nbpoints(isub),isub)=X(ind_sel(ind_ind_sel)); |
---|
| 516 | Y_tps(1:nbpoints(isub),isub)=Y(ind_sel(ind_ind_sel)); |
---|
| 517 | U_tps(1:nbpoints(isub)+3,isub)=U_tps_sub; |
---|
| 518 | V_tps(1:nbpoints(isub)+3,isub)=V_tps_sub; |
---|
| 519 | nb_select(ind_sel(ind_ind_sel))=nb_select(ind_sel(ind_ind_sel))+1; |
---|
| 520 | display('good2') |
---|
| 521 | break |
---|
| 522 | end |
---|
| 523 | end |
---|
[242] | 524 | end |
---|
| 525 | end |
---|
[246] | 526 | ind_empty=find(test_empty); |
---|
| 527 | %remove empty subdomains |
---|
| 528 | if ~isempty(ind_empty) |
---|
| 529 | SubRangx(ind_empty,:)=[]; |
---|
| 530 | SubRangy(ind_empty,:)=[]; |
---|
| 531 | X_tps(:,ind_empty)=[]; |
---|
| 532 | Y_tps(:,ind_empty)=[]; |
---|
| 533 | U_tps(:,ind_empty)=[]; |
---|
| 534 | V_tps(:,ind_empty)=[]; |
---|
| 535 | end |
---|
| 536 | nb_select(nb_select==0)=1;%ones(size(find(nb_select==0))); |
---|
| 537 | U_smooth=U_smooth./nb_select; |
---|
| 538 | V_smooth=V_smooth./nb_select; |
---|
[242] | 539 | |
---|
| 540 | |
---|
| 541 | |
---|
[246] | 542 | |
---|
| 543 | |
---|
| 544 | |
---|
[242] | 545 | % U_patch = EM * spline_coeff; |
---|
| 546 | % U_patch=reshape(U_patch,npy,npx); |
---|
| 547 | % PM = [ones(size(dsites,1),1) dsites]; |
---|
| 548 | % EM = [IM_sites PM]; |
---|
| 549 | % U(test_false)=[]; |
---|
| 550 | % U_nodes=EM * spline_coeff; |
---|
| 551 | |
---|
| 552 | %exact = testfunctions(epoints); |
---|
| 553 | %maxerr = norm(Pf-exact,inf); |
---|
| 554 | % PlotSurf(xe,ye,Pf,neval,exact,maxerr,[160,20]); |
---|
| 555 | % PlotError2D(xe,ye,Pf,exact,maxerr,neval,[160,20]); |
---|
| 556 | |
---|
| 557 | |
---|
| 558 | |
---|
| 559 | % DM = DistanceMatrix(dsites,ctrs) |
---|
| 560 | % Forms the distance matrix of two sets of points in R^s, |
---|
| 561 | % i.e., DM(i,j) = || datasite_i - center_j ||_2. |
---|
| 562 | % Input |
---|
| 563 | % dsites: Mxs matrix representing a set of M data sites in R^s |
---|
| 564 | % (i.e., each row contains one s-dimensional point) |
---|
| 565 | % ctrs: Nxs matrix representing a set of N centers in R^s |
---|
| 566 | % (one center per row) |
---|
| 567 | % Output |
---|
| 568 | |
---|
| 569 | |
---|
| 570 | |
---|
| 571 | |
---|