source: trunk/src/civ_uvmat.m @ 292

Last change on this file since 292 was 273, checked in by sommeria, 12 years ago

GUI civ renovated with panels

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