source: trunk/src/civ_matlab.m @ 589

Last change on this file since 589 was 582, checked in by sommeria, 11 years ago

various bugs corrected

File size: 34.3 KB
Line 
1%'civ_matlab': main PIV function, calleed by the GUI civ
2% --- call the sub-functions:
3%   civ: PIV function itself
4%   fix: removes false vectors after detection by various criteria
5%   filter_tps: make interpolation-smoothing
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%     .Civ1: for civ1
17%     .Fix1:
18%     .Patch1:
19%     .Civ2: for civ2
20%     .Fix2:
21%     .Patch2:
22% ncfile: name of a netcdf file to be created for the result (extension .nc)
23%
24%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
25%  Copyright  2011, LEGI / CNRS-UJF-INPG, joel.sommeria@legi.grenoble-inp.fr.
26%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
27%     This is part of the toolbox UVMAT.
28%
29%     UVMAT is free software; you can redistribute it and/or modify
30%     it under the terms of the GNU General Public License as published by
31%     the Free Software Foundation; either version 2 of the License, or
32%     (at your option) any later version.
33%
34%     UVMAT is distributed in the hope that it will be useful,
35%     but WITHOUT ANY WARRANTY; without even the implied warranty of
36%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
37%     GNU General Public License (open UVMAT/COPYING.txt) for more details.
38%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
39
40function [Data,errormsg,result_conv]= civ_matlab(Param,ncfile)
41errormsg='';
42Data.ListGlobalAttribute={'Conventions','Program','CivStage'};
43Data.Conventions='uvmat/civdata';% states the conventions used for the description of field variables and attributes
44Data.Program='civ_matlab';
45Data.CivStage=0;%default
46ListVarCiv1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'}; %variables to read
47ListVarFix1={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F','Civ1_FF'};
48mask='';
49maskname='';%default
50check_civx=0;%default
51check_civ1=0;%default
52check_patch1=0;%default
53
54% case of input Param set by an xml file (batch mode)
55if ischar(Param)
56    Param=xml2struct(Param); %if Param is the name of an xml file, read this file as a Matlab structure
57end
58
59%% Civ1
60if isfield (Param,'Civ1')
61    par_civ1=Param.Civ1;
62    if isfield(par_civ1,'reverse_pair')% A REVOIR
63        if par_civ1.reverse_pair
64            if ischar(par_civ1.ImageB)
65                temp=par_civ1.ImageA;
66                par_civ1.ImageA=imread(par_civ1.ImageB);
67            end
68            if ischar(temp)
69                par_civ1.ImageB=imread(temp);
70            end
71        end
72    else
73        if isfield(Param.Civ1,'ImageA') && ischar(Param.Civ1.ImageA)
74             Param.Civ1.ImageA=regexprep(Param.Civ1.ImageA,'''','\');
75            [par_civ1.ImageA,VideoObject] = read_image(Param.Civ1.ImageA,par_civ1.FileTypeA,[],par_civ1.FrameIndexA);
76        end
77        if isfield(Param.Civ1,'ImageB')&& ischar(Param.Civ1.ImageB)
78             Param.Civ1.ImageB=regexprep(Param.Civ1.ImageB,'''','\');
79             if strcmp(Param.Civ1.ImageA,Param.Civ1.ImageB)% use the same movie object
80                 [par_civ1.ImageB,VideoObject] = read_image(Param.Civ1.ImageB,par_civ1.FileTypeB,VideoObject,par_civ1.FrameIndexB);
81             else
82            [par_civ1.ImageB,VideoObject] = read_image(Param.Civ1.ImageB,par_civ1.FileTypeB,par_civ1.ImageB,par_civ1.FrameIndexB);
83             end
84        end
85    end
86   
87    % caluclate velocity data (y and v in indices, reverse to y component)
88    [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
89    if ~isempty(errormsg)
90        return
91    end
92    list_param=(fieldnames(Param.Civ1))';
93    Civ1_param=list_param;%default
94    %set the values of all the global attributes in list_param
95    for ilist=1:length(list_param)
96        Civ1_param{ilist}=['Civ1_' list_param{ilist}];
97        Data.(['Civ1_' list_param{ilist}])=Param.Civ1.(list_param{ilist});
98    end
99    Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param];% {'Civ1_Time','Civ1_Dt'}];
100    Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_F','Civ1_C'};%  cell array containing the names of the fields to record
101    Data.VarDimName={'nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1','nb_vec_1'};
102    Data.VarAttribute{1}.Role='coord_x';
103    Data.VarAttribute{2}.Role='coord_y';
104    Data.VarAttribute{3}.Role='vector_x';
105    Data.VarAttribute{4}.Role='vector_y';
106    Data.VarAttribute{5}.Role='warnflag';
107    Data.Civ1_X=reshape(xtable,[],1);
108    Data.Civ1_Y=reshape(Param.Civ1.ImageHeight-ytable+1,[],1);
109    Data.Civ1_U=reshape(utable,[],1);
110    Data.Civ1_V=reshape(-vtable,[],1);
111    Data.Civ1_C=reshape(ctable,[],1);
112    Data.Civ1_F=reshape(F,[],1);
113    Data.CivStage=1;
114else
115    if exist('ncfile','var')
116        CivFile=ncfile;
117    elseif isfield(Param.Patch1,'CivFile')
118        CivFile=Param.Patch1.CivFile;
119    end
120    Data=nc2struct(CivFile,'ListGlobalAttribute','absolut_time_T0'); %look for the constant 'absolut_time_T0' to detect old civx data format
121    if isfield(Data,'Txt')
122        errormsg=Data.Txt;
123        return
124    end
125    if ~isempty(Data.absolut_time_T0')%read civx file
126        check_civx=1;% test for old civx data format
127        [Data,vardetect,ichoice]=nc2struct(CivFile);%read the variables in the netcdf file
128    else
129        Data=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
130    end
131end
132
133%% Fix1
134if isfield (Param,'Fix1')
135    ListFixParam=fieldnames(Param.Fix1);
136    for ilist=1:length(ListFixParam)
137        ParamName=ListFixParam{ilist};
138        ListName=['Fix1_' ParamName];
139        eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
140        eval(['Data.' ListName '=Param.Fix1.' ParamName ';'])
141    end
142    if check_civx
143        if ~isfield(Data,'fix')
144            Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix'];
145            Data.fix=1;
146            Data.ListVarName=[Data.ListVarName {'vec_FixFlag'}];
147            Data.VarDimName=[Data.VarDimName {'nb_vectors'}];
148        end
149        Data.vec_FixFlag=fix(Param.Fix1,Data.vec_F,Data.vec_C,Data.vec_U,Data.vec_V,Data.vec_X,Data.vec_Y);
150    else
151        Data.ListVarName=[Data.ListVarName {'Civ1_FF'}];
152        Data.VarDimName=[Data.VarDimName {'nb_vec_1'}];
153        nbvar=length(Data.ListVarName);
154        Data.VarAttribute{nbvar}.Role='errorflag';
155        Data.Civ1_FF=fix(Param.Fix1,Data.Civ1_F,Data.Civ1_C,Data.Civ1_U,Data.Civ1_V);
156        Data.CivStage=2;
157    end
158end
159%% Patch1
160if isfield (Param,'Patch1')
161    if check_civx
162        errormsg='Civ Matlab input needed for patch';
163        return
164    end
165   
166    Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch1_Rho','Patch1_Threshold','Patch1_SubDomain'}];
167    Data.Patch1_FieldSmooth=Param.Patch1.FieldSmooth;
168    Data.Patch1_MaxDiff=Param.Patch1.MaxDiff;
169    Data.Patch1_SubDomainSize=Param.Patch1.SubDomainSize;
170    nbvar=length(Data.ListVarName);
171    Data.ListVarName=[Data.ListVarName {'Civ1_U_smooth','Civ1_V_smooth','Civ1_SubRange','Civ1_NbCentres','Civ1_Coord_tps','Civ1_U_tps','Civ1_V_tps'}];
172    Data.VarDimName=[Data.VarDimName {'nb_vec_1','nb_vec_1',{'nb_coord','nb_bounds','nb_subdomain_1'},'nb_subdomain_1',...
173        {'nb_tps_1','nb_coord','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'},{'nb_tps_1','nb_subdomain_1'}}];
174    Data.VarAttribute{nbvar+1}.Role='vector_x';
175    Data.VarAttribute{nbvar+2}.Role='vector_y';
176    Data.VarAttribute{nbvar+5}.Role='coord_tps';
177    Data.VarAttribute{nbvar+6}.Role='vector_x';
178    Data.VarAttribute{nbvar+7}.Role='vector_y';
179    Data.Civ1_U_smooth=zeros(size(Data.Civ1_X));
180    Data.Civ1_V_smooth=zeros(size(Data.Civ1_X));
181    if isfield(Data,'Civ1_FF')
182        ind_good=find(Data.Civ1_FF==0);
183    else
184        ind_good=1:numel(Data.Civ1_X);
185    end
186    [Data.Civ1_SubRange,Data.Civ1_NbCentres,Data.Civ1_Coord_tps,Data.Civ1_U_tps,Data.Civ1_V_tps,tild,Ures, Vres,tild,FFres]=...
187        filter_tps([Data.Civ1_X(ind_good) Data.Civ1_Y(ind_good)],Data.Civ1_U(ind_good),Data.Civ1_V(ind_good),[],Data.Patch1_SubDomainSize,Data.Patch1_FieldSmooth,Data.Patch1_MaxDiff);
188    Data.Civ1_U_smooth(ind_good)=Ures;
189    Data.Civ1_V_smooth(ind_good)=Vres;
190    Data.Civ1_FF(ind_good)=FFres;
191    Data.CivStage=3;
192end
193
194%% Civ2
195if isfield (Param,'Civ2')
196    par_civ2=Param.Civ2;
197    par_civ2.ImageA=[];
198    par_civ2.ImageB=[];
199    if isfield(Param.Civ2,'ImageA') && isfield(Param.Civ2,'ImageB')
200        Param.Civ2.ImageA=regexprep(Param.Civ2.ImageA,'''','\');
201        Param.Civ2.ImageB=regexprep(Param.Civ2.ImageB,'''','\');
202        if isfield (Param,'Civ1')
203            Param.Civ2.ImageA=regexprep(Param.Civ2.ImageA,'''','\');
204            if strcmp(Param.Civ1.ImageA,Param.Civ2.ImageA)
205                if isequal(Param.Civ1.FrameIndexA,Param.Civ2.FrameIndexA)
206                    par_civ2.ImageA=par_civ1.ImageA;
207                else % read another frame of the same movie object
208                    par_civ2.ImageA = read_image(Param.Civ2.ImageA,Param.Civ2.FileTypeA,VideoObject,Param.Civ2.FrameIndexA);
209                end
210            end
211            Param.Civ2.ImageB=regexprep(Param.Civ2.ImageB,'''','\');
212            if strcmp(Param.Civ1.ImageB,Param.Civ2.ImageB)
213                if isequal(Param.Civ1.FrameIndexB,Param.Civ2.FrameIndexB)
214                    par_civ2.ImageB=par_civ1.ImageB;
215                else % read another frame of the same movie object
216                    par_civ2.ImageB = read_image(Param.Civ2.ImageB,Param.Civ2.FileTypeB,VideoObject,Param.Civ2.FrameIndexB);
217                end
218            end
219        end
220        if isempty(par_civ2.ImageA) && isfield(Param.Civ2,'ImageA')
221            [par_civ2.ImageA,VideoObject] = read_image(Param.Civ2.ImageA,Param.Civ2.FileTypeA,[],Param.Civ2.FrameIndexA);
222        end
223        if isempty(par_civ2.ImageB)&& isfield(Param.Civ2,'ImageB')
224            if strcmp(Param.Civ2.ImageA,Param.Civ2.ImageB)
225                par_civ2.ImageB = read_image(Param.Civ2.ImageB,Param.Civ2.FileTypeB,VideoObject,Param.Civ2.FrameIndexB);
226            else
227                par_civ2.ImageB = read_image(Param.Civ2.ImageB,Param.Civ2.FileTypeB,[],Param.Civ2.FrameIndexB);
228            end
229        end
230    end
231    %     ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
232    %     iby2=ceil(par_civ2.CorrBoxSize(2)/2);
233    %     isx2=ibx2+4;% search ara +-4 pixels around the guess
234    %     isy2=iby2+4;
235    %     % shift from par_civ2.filename_nc1
236    %     % shiftx=velocity interpolated at position
237    %     miniy=max(1+isy2,1+iby2);
238    %     minix=max(1+isx2,1+ibx2);
239    %     maxiy=min(size(par_civ2.ImageA,1)-isy2,size(par_civ2.ImageA,1)-iby2);
240    %     maxix=min(size(par_civ2.ImageA,2)-isx2,size(par_civ2.ImageA,2)-ibx2);
241    %     [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
242    %     GridX=reshape(GridX,[],1);
243    %     GridY=reshape(GridY,[],1);
244    if isfield(par_civ2,'Grid')% grid points set as input file
245        if ischar(par_civ2.Grid)%read the grid file if the input is a file name
246            par_civ2.Grid=dlmread(par_civ2.Grid);
247            par_civ2.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
248        end
249    else% automatic grid
250        minix=floor(par_civ2.Dx/2)-0.5;
251        maxix=minix+par_civ2.Dx*floor((par_civ2.ImageWidth-1)/par_civ2.Dx);
252        miniy=floor(par_civ2.Dy/2)-0.5;
253        maxiy=minix+par_civ2.Dy*floor((par_civ2.ImageHeight-1)/par_civ2.Dy);
254        [GridX,GridY]=meshgrid(minix:par_civ2.Dx:maxix,miniy:par_civ2.Dy:maxiy);
255        par_civ2.Grid(:,1)=reshape(GridX,[],1);
256        par_civ2.Grid(:,2)=reshape(GridY,[],1);
257    end
258    Shiftx=zeros(size(par_civ2.Grid,1),1);% shift expected from civ1 data
259    Shifty=zeros(size(par_civ2.Grid,1),1);
260    nbval=zeros(size(par_civ2.Grid,1),1);
261    if par_civ2.CheckDeformation
262        DUDX=zeros(size(par_civ2.Grid,1),1);
263        DUDY=zeros(size(par_civ2.Grid,1),1);
264        DVDX=zeros(size(par_civ2.Grid,1),1);
265        DVDY=zeros(size(par_civ2.Grid,1),1);
266    end
267    NbSubDomain=size(Data.Civ1_SubRange,3);
268    % get the guess from patch1
269    for isub=1:NbSubDomain
270        nbvec_sub=Data.Civ1_NbCentres(isub);
271        ind_sel=find(GridX>=Data.Civ1_SubRange(1,1,isub) & GridX<=Data.Civ1_SubRange(1,2,isub) & GridY>=Data.Civ1_SubRange(2,1,isub) & GridY<=Data.Civ1_SubRange(2,2,isub));
272        epoints = [GridX(ind_sel) GridY(ind_sel)];% coordinates of interpolation sites
273        ctrs=Data.Civ1_Coord_tps(1:nbvec_sub,:,isub) ;%(=initial points) ctrs
274        nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
275        EM = tps_eval(epoints,ctrs);
276        Shiftx(ind_sel)=Shiftx(ind_sel)+EM*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
277        Shifty(ind_sel)=Shifty(ind_sel)+EM*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
278        if par_civ2.CheckDeformation
279            [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
280            DUDX(ind_sel)=DUDX(ind_sel)+EMDX*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
281            DUDY(ind_sel)=DUDY(ind_sel)+EMDY*Data.Civ1_U_tps(1:nbvec_sub+3,isub);
282            DVDX(ind_sel)=DVDX(ind_sel)+EMDX*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
283            DVDY(ind_sel)=DVDY(ind_sel)+EMDY*Data.Civ1_V_tps(1:nbvec_sub+3,isub);
284        end
285    end
286    mask='';
287    if par_civ2.CheckMask&&~isempty(par_civ2.Mask)&& ~strcmp(maskname,par_civ2.Mask)% mask exist, not already read in civ1
288        mask=imread(par_civ2.Mask);
289    end
290    ibx2=ceil(par_civ2.CorrBoxSize(1)/2);
291    iby2=ceil(par_civ2.CorrBoxSize(2)/2);
292    %     isx2=ibx2+4;% search ara +-4 pixels around the guess
293    %     isy2=iby2+4;
294    par_civ2.SearchBoxSize(1)=2*ibx2+9;% search ara +-4 pixels around the guess
295    par_civ2.SearchBoxSize(2)=2*iby2+9;
296    %par_civ2.SearchBoxSize(1)=2*isx2+1;
297    %par_civ2.SearchBoxSize(2)=2*isy2+1;
298    par_civ2.SearchBoxShift=[Shiftx(nbval>=1)./nbval(nbval>=1) Shifty(nbval>=1)./nbval(nbval>=1)];
299    par_civ2.Grid=[GridX(nbval>=1)-par_civ2.SearchBoxShift(:,1)/2 GridY(nbval>=1)-par_civ2.SearchBoxShift(:,2)/2];% grid taken at the extrapolated origin of the displacement vectors
300    if par_civ2.CheckDeformation
301        par_civ2.DUDX=DUDX./nbval;
302        par_civ2.DUDY=DUDY./nbval;
303        par_civ2.DVDX=DVDX./nbval;
304        par_civ2.DVDY=DVDY./nbval;
305    end
306    % caluclate velocity data (y and v in indices, reverse to y component)
307    [xtable ytable utable vtable ctable F] = civ (par_civ2);
308    list_param=(fieldnames(Param.Civ2))';
309    list_remove={'pxcmx','pxcmy','npx','npy','gridflag','maskflag','term_a','term_b','T0'};
310    for ilist=1:length(list_remove)
311        index=strcmp(list_remove{ilist},list_param);
312        if ~isempty(find(index,1))
313            list_param(index)=[];
314        end
315    end
316    for ilist=1:length(list_param)
317        Civ2_param{ilist}=['Civ2_' list_param{ilist}];
318        eval(['Data.Civ2_' list_param{ilist} '=Param.Civ2.' list_param{ilist} ';'])
319    end
320    if isfield(Data,'Civ2_gridname') && strcmp(Data.Civ1_gridname(1:6),'noFile')
321        Data.Civ1_gridname='';
322    end
323    if isfield(Data,'Civ2_maskname') && strcmp(Data.Civ1_maskname(1:6),'noFile')
324        Data.Civ2_maskname='';
325    end
326    Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ2_param {'Civ2_Time','Civ2_Dt'}];
327    nbvar=numel(Data.ListVarName);
328    Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];%  cell array containing the names of the fields to record
329    Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2','nb_vec_2'}];
330    Data.VarAttribute{nbvar+1}.Role='coord_x';
331    Data.VarAttribute{nbvar+2}.Role='coord_y';
332    Data.VarAttribute{nbvar+3}.Role='vector_x';
333    Data.VarAttribute{nbvar+4}.Role='vector_y';
334    Data.VarAttribute{nbvar+5}.Role='warnflag';
335    Data.Civ2_X=reshape(xtable,[],1);
336    Data.Civ2_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
337    Data.Civ2_U=reshape(utable,[],1);
338    Data.Civ2_V=reshape(-vtable,[],1);
339    Data.Civ2_C=reshape(ctable,[],1);
340    Data.Civ2_F=reshape(F,[],1);
341    Data.CivStage=Data.CivStage+1;
342end
343
344%% Fix2
345if isfield (Param,'Fix2')
346    ListFixParam=fieldnames(Param.Fix2);
347    for ilist=1:length(ListFixParam)
348        ParamName=ListFixParam{ilist};
349        ListName=['Fix2_' ParamName];
350        eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
351        eval(['Data.' ListName '=Param.Fix2.' ParamName ';'])
352    end
353    if check_civx
354        if ~isfield(Data,'fix2')
355            Data.ListGlobalAttribute=[Data.ListGlobalAttribute 'fix2'];
356            Data.fix2=1;
357            Data.ListVarName=[Data.ListVarName {'vec2_FixFlag'}];
358            Data.VarDimName=[Data.VarDimName {'nb_vectors2'}];
359        end
360        Data.vec_FixFlag=fix(Param.Fix2,Data.vec2_F,Data.vec2_C,Data.vec2_U,Data.vec2_V,Data.vec2_X,Data.vec2_Y);
361    else
362        Data.ListVarName=[Data.ListVarName {'Civ2_FF'}];
363        Data.VarDimName=[Data.VarDimName {'nb_vec_2'}];
364        nbvar=length(Data.ListVarName);
365        Data.VarAttribute{nbvar}.Role='errorflag';
366        Data.Civ2_FF=fix(Param.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V);
367        Data.CivStage=Data.CivStage+1;
368    end
369   
370end
371
372%% Patch2
373if isfield (Param,'Patch2')
374    Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch2_Rho','Patch2_Threshold','Patch2_SubDomain'}];
375    Data.Patch2_FieldSmooth=Param.Patch2.FieldSmooth;
376    Data.Patch2_MaxDiff=Param.Patch2.MaxDiff;
377    Data.Patch2_SubDomainSize=Param.Patch2.SubDomainSize;
378    nbvar=length(Data.ListVarName);
379    Data.ListVarName=[Data.ListVarName {'Civ2_U_smooth','Civ2_V_smooth','Civ2_SubRange','Civ2_NbCentres','Civ2_Coord_tps','Civ2_U_tps','Civ2_V_tps'}];
380    Data.VarDimName=[Data.VarDimName {'nb_vec_2','nb_vec_2',{'nb_coord','nb_bounds','nb_subdomain_2'},{'nb_subdomain_2'},...
381        {'nb_tps_2','nb_coord','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'},{'nb_tps_2','nb_subdomain_2'}}];
382   
383    Data.VarAttribute{nbvar+1}.Role='vector_x';
384    Data.VarAttribute{nbvar+2}.Role='vector_y';
385    Data.VarAttribute{nbvar+5}.Role='coord_tps';
386    Data.VarAttribute{nbvar+6}.Role='vector_x';
387    Data.VarAttribute{nbvar+7}.Role='vector_y';
388    Data.Civ2_U_smooth=zeros(size(Data.Civ2_X));
389    Data.Civ2_V_smooth=zeros(size(Data.Civ2_X));
390    if isfield(Data,'Civ2_FF')
391        ind_good=find(Data.Civ2_FF==0);
392    else
393        ind_good=1:numel(Data.Civ2_X);
394    end
395    [Data.Civ2_SubRange,Data.Civ2_NbCentres,Data.Civ2_Coord_tps,Data.Civ2_U_tps,Data.Civ2_V_tps,tild,Ures, Vres,tild,FFres]=...
396        filter_tps([Data.Civ2_X(ind_good) Data.Civ2_Y(ind_good)],Data.Civ2_U(ind_good),Data.Civ2_V(ind_good),[],Data.Patch2_SubDomainSize,Data.Patch2_FieldSmooth,Data.Patch2_MaxDiff);
397    Data.Civ2_U_smooth(ind_good)=Ures;
398    Data.Civ2_V_smooth(ind_good)=Vres;
399    Data.Civ2_FF(ind_good)=FFres;
400    Data.CivStage=Data.CivStage+1;
401end
402
403%% write result in a netcdf file if requested
404if exist('ncfile','var')
405    errormsg=struct2nc(ncfile,Data);
406    if isempty(errormsg)
407        disp([ncfile ' written'])
408    else
409        disp(errormsg)
410    end
411end
412
413% 'civ': function piv.m adapted from PIVlab http://pivlab.blogspot.com/
414%--------------------------------------------------------------------------
415% function [xtable ytable utable vtable typevector] = civ (image1,image2,ibx,iby step, subpixfinder, mask, roi)
416%
417% OUTPUT:
418% xtable: set of x coordinates
419% ytable: set of y coordiantes
420% utable: set of u displacements (along x)
421% vtable: set of v displacements (along y)
422% ctable: max image correlation for each vector
423% typevector: set of flags, =1 for good, =0 for NaN vectors
424%
425%INPUT:
426% par_civ: structure of input parameters, with fields:
427%  .CorrBoxSize
428%  .SearchBoxSize
429%  .SearchBoxShift
430%  .ImageHeight
431%  .ImageWidth
432%  .Dx, Dy
433%  .Grid
434%  .Mask
435%  .MinIma
436%  .MaxIma
437%  .image1:first image (matrix)
438% image2: second image (matrix)
439% ibx2,iby2: half size of the correlation box along x and y, in px (size=(2*iby2+1,2*ibx2+1)
440% isx2,isy2: half size of the search box along x and y, in px (size=(2*isy2+1,2*isx2+1)
441% shiftx, shifty: shift of the search box (in pixel index, yshift reversed)
442% step: mesh of the measurement points (in px)
443% subpixfinder=1 or 2 controls the curve fitting of the image correlation
444% mask: =[] for no mask
445% roi: 4 element vector defining a region of interest: x position, y position, width, height, (in image indices), for the whole image, roi=[];
446function [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ)
447%this funtion performs the DCC PIV analysis. Recent window-deformation
448%methods perform better and will maybe be implemented in the future.
449
450%% prepare measurement grid
451if isfield(par_civ,'Grid')% grid points set as input
452    if ischar(par_civ.Grid)%read the drid file if the input is a file name
453        par_civ.Grid=dlmread(par_civ.Grid);
454        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
455    end
456else% automatic grid
457    minix=floor(par_civ.Dx/2)-0.5;
458    maxix=minix+par_civ.Dx*floor((par_civ.ImageWidth-1)/par_civ.Dx);
459    miniy=floor(par_civ.Dy/2)-0.5;
460    maxiy=minix+par_civ.Dy*floor((par_civ.ImageHeight-1)/par_civ.Dy);
461    [GridX,GridY]=meshgrid(minix:par_civ.Dx:maxix,miniy:par_civ.Dy:maxiy);
462    par_civ.Grid(:,1)=reshape(GridX,[],1);
463    par_civ.Grid(:,2)=reshape(GridY,[],1);
464end
465nbvec=size(par_civ.Grid,1);
466
467%% prepare correlation and search boxes
468ibx2=ceil(par_civ.CorrBoxSize(1)/2);
469iby2=ceil(par_civ.CorrBoxSize(2)/2);
470isx2=ceil(par_civ.SearchBoxSize(1)/2);
471isy2=ceil(par_civ.SearchBoxSize(2)/2);
472shiftx=round(par_civ.SearchBoxShift(:,1));
473shifty=-round(par_civ.SearchBoxShift(:,2));% sign minus because image j index increases when y decreases
474if numel(shiftx)==1% case of a unique shift for the whole field( civ1)
475    shiftx=shiftx*ones(nbvec,1);
476    shifty=shifty*ones(nbvec,1);
477end
478
479%% Default output
480xtable=par_civ.Grid(:,1);
481ytable=par_civ.Grid(:,2);
482utable=zeros(nbvec,1);
483vtable=zeros(nbvec,1);
484ctable=zeros(nbvec,1);
485F=zeros(nbvec,1);
486result_conv=[];
487errormsg='';
488
489%% prepare mask
490if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
491    if strcmp(par_civ.Mask,'all')
492        return    % get the grid only, no civ calculation
493    elseif ischar(par_civ.Mask)
494        par_civ.Mask=imread(par_civ.Mask);
495    end
496end
497check_MinIma=isfield(par_civ,'MinIma');% test for image luminosity threshold
498check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
499
500% %% prepare images
501% if isfield(par_civ,'reverse_pair')
502%     if par_civ.reverse_pair
503%         if ischar(par_civ.ImageB)
504%             temp=par_civ.ImageA;
505%             par_civ.ImageA=imread(par_civ.ImageB);
506%         end
507%         if ischar(temp)
508%             par_civ.ImageB=imread(temp);
509%         end
510%     end
511% else
512%     if ischar(par_civ.ImageA)
513%         par_civ.ImageA=imread(par_civ.ImageA);
514%     end
515%     if ischar(par_civ.ImageB)
516%         par_civ.ImageB=imread(par_civ.ImageB);
517%     end
518% end
519par_civ.ImageA=sum(double(par_civ.ImageA),3);%sum over rgb component for color images
520par_civ.ImageB=sum(double(par_civ.ImageB),3);
521[npy_ima npx_ima]=size(par_civ.ImageA);
522if ~isequal(size(par_civ.ImageB),[npy_ima npx_ima])
523    errormsg='image pair with unequal size';
524    return
525end
526
527%% Apply mask
528    % Convention for mask IDEAS TO IMPLEMENT ?
529    % mask >200 : velocity calculated
530    %  200 >=mask>150;velocity not calculated, interpolation allowed (bad spots)
531    % 150>=mask >100: velocity not calculated, nor interpolated
532    %  100>=mask> 20: velocity not calculated, impermeable (no flux through mask boundaries)
533    %  20>=mask: velocity=0
534checkmask=0;
535MinA=min(min(par_civ.ImageA));
536MinB=min(min(par_civ.ImageB));
537if isfield(par_civ,'Mask') && ~isempty(par_civ.Mask)
538   checkmask=1;
539   if ~isequal(size(par_civ.Mask),[npy_ima npx_ima])
540        errormsg='mask must be an image with the same size as the images';
541        return
542   end
543  %  check_noflux=(par_civ.Mask<100) ;%TODO: to implement
544    check_undefined=(par_civ.Mask<200 & par_civ.Mask>=20 );
545    par_civ.ImageA(check_undefined)=MinA;% put image A to zero (i.e. the min image value) in the undefined  area
546    par_civ.ImageB(check_undefined)=MinB;% put image B to zero (i.e. the min image value) in the undefined  area
547end
548
549%% compute image correlations: MAINLOOP on velocity vectors
550corrmax=0;
551sum_square=1;% default
552mesh=1;% default
553CheckDecimal=isfield(par_civ,'CheckDecimal')&& par_civ.CheckDecimal==1;
554if CheckDecimal
555    mesh=0.2;%mesh in pixels for subpixel image interpolation
556    CheckDeformation=isfield(par_civ,'CheckDeformation')&& par_civ.CheckDeformation==1;
557end
558% vector=[0 0];%default
559for ivec=1:nbvec
560    iref=round(par_civ.Grid(ivec,1)+0.5);% xindex on the image A for the middle of the correlation box
561    jref=round(par_civ.ImageHeight-par_civ.Grid(ivec,2)+0.5);% yindex on the image B for the middle of the correlation box
562    if ~(checkmask && par_civ.Mask(jref,iref)<=20) %velocity not set to zero by the black mask
563%         if jref-iby2<1 || jref+iby2>par_civ.ImageHeight|| iref-ibx2<1 || iref+ibx2>par_civ.ImageWidth||...
564%               jref+shifty(ivec)-isy2<1||jref+shifty(ivec)+isy2>par_civ.ImageHeight|| iref+shiftx(ivec)-isx2<1 || iref+shiftx(ivec)+isx2>par_civ.ImageWidth  % we are outside the image
565%             F(ivec)=3;
566%         else
567            F(ivec)=0;
568            subrange1_x=iref-ibx2:iref+ibx2;% x indices defining the first subimage
569            subrange1_y=jref-iby2:jref+iby2;% y indices defining the first subimage
570            subrange2_x=iref+shiftx(ivec)-isx2:iref+shiftx(ivec)+isx2;%x indices defining the second subimage
571            subrange2_y=jref+shifty(ivec)-isy2:jref+shifty(ivec)+isy2;%y indices defining the second subimage
572            image1_crop=MinA*ones(numel(subrange1_y),numel(subrange1_x));% default value=min of image A
573            image2_crop=MinA*ones(numel(subrange2_y),numel(subrange2_x));% default value=min of image A
574            check1_x=subrange1_x>=1 & subrange1_x<=par_civ.ImageWidth;% check which points in the subimage 1 are contained in the initial image 1
575            check1_y=subrange1_y>=1 & subrange1_y<=par_civ.ImageHeight;
576            check2_x=subrange2_x>=1 & subrange2_x<=par_civ.ImageWidth;% check which points in the subimage 2 are contained in the initial image 2
577            check2_y=subrange2_y>=1 & subrange2_y<=par_civ.ImageHeight;
578            image1_crop(check1_y,check1_x)=par_civ.ImageA(subrange1_y(check1_y),subrange1_x(check1_x));%extract a subimage (correlation box) from image A
579            image2_crop(check2_y,check2_x)=par_civ.ImageB(subrange2_y(check2_y),subrange2_x(check2_x));%extract a larger subimage (search box) from image B
580            image1_mean=mean(mean(image1_crop));
581            image2_mean=mean(mean(image2_crop));
582            %threshold on image minimum
583            if check_MinIma && (image1_mean < par_civ.MinIma || image2_mean < par_civ.MinIma)
584                F(ivec)=3;
585            end
586            %threshold on image maximum
587            if check_MaxIma && (image1_mean > par_civ.MaxIma || image2_mean > par_civ.MaxIma)
588                F(ivec)=3;
589            end
590%         end     
591        if F(ivec)~=3
592            image1_crop=image1_crop-image1_mean;%substract the mean
593            image2_crop=image2_crop-image2_mean;
594            if CheckDecimal
595                xi=(1:mesh:size(image1_crop,2));
596                yi=(1:mesh:size(image1_crop,1))';
597                if CheckDeformation
598                    [XI,YI]=meshgrid(xi-ceil(size(image1_crop,2)/2),yi-ceil(size(image1_crop,1)/2));
599                    XIant=XI-par_civ.DUDX(ivec)*XI-par_civ.DUDY(ivec)*YI+ceil(size(image1_crop,2)/2);
600                    YIant=YI-par_civ.DVDX(ivec)*XI-par_civ.DVDY(ivec)*YI+ceil(size(image1_crop,1)/2);
601                    image1_crop=interp2(image1_crop,XIant,YIant);
602                else
603                    image1_crop=interp2(image1_crop,xi,yi);
604                end
605                xi=(1:mesh:size(image2_crop,2));
606                yi=(1:mesh:size(image2_crop,1))';
607                image2_crop=interp2(image2_crop,xi,yi);
608            end
609            sum_square=sum(sum(image1_crop.*image1_crop));
610            %reference: Oliver Pust, PIV: Direct Cross-Correlation
611            result_conv= conv2(image2_crop,flipdim(flipdim(image1_crop,2),1),'valid');
612            corrmax= max(max(result_conv));
613            result_conv=(result_conv/corrmax)*255; %normalize, peak=always 255
614            %Find the correlation max, at 255
615            [y,x] = find(result_conv==255,1);
616            if ~isempty(y) && ~isempty(x)
617                try
618                    if par_civ.CorrSmooth==1
619                        [vector,F(ivec)] = SUBPIXGAUSS (result_conv,x,y);
620                    elseif par_civ.CorrSmooth==2
621                        [vector,F(ivec)] = SUBPIX2DGAUSS (result_conv,x,y);
622                    end
623                    utable(ivec)=vector(1)*mesh+shiftx(ivec);
624                    vtable(ivec)=vector(2)*mesh+shifty(ivec);                 
625                    xtable(ivec)=iref+utable(ivec)/2-0.5;% convec flow (velocity taken at the point middle from imgae 1 and 2)
626                    ytable(ivec)=jref+vtable(ivec)/2-0.5;% and position of pixel 1=0.5 (convention for image coordinates=0 at the edge)
627                    iref=round(xtable(ivec));% image index for the middle of the vector
628                    jref=round(ytable(ivec));
629                    if checkmask && par_civ.Mask(jref,iref)<200 && par_civ.Mask(jref,iref)>=100
630                        utable(ivec)=0;
631                        vtable(ivec)=0;
632                        F(ivec)=3;
633                    end
634                    ctable(ivec)=corrmax/sum_square;% correlation value
635                catch ME
636                    F(ivec)=3;
637                end
638            else
639                F(ivec)=3;
640            end
641        end
642    end
643   
644    %Create the vector matrix x, y, u, v
645end
646result_conv=result_conv*corrmax/(255*sum_square);% keep the last correlation matrix for output
647
648%------------------------------------------------------------------------
649% --- Find the maximum of the correlation function after interpolation
650function [vector,F] = SUBPIXGAUSS (result_conv,x,y)
651%------------------------------------------------------------------------
652vector=[0 0]; %default
653F=0;
654[npy,npx]=size(result_conv);
655
656% if (x <= (size(result_conv,1)-1)) && (y <= (size(result_conv,1)-1)) && (x >= 1) && (y >= 1)
657    %the following 8 lines are copyright (c) 1998, Uri Shavit, Roi Gurka, Alex Liberzon, Technion ï¿œ Israel Institute of Technology
658    %http://urapiv.wordpress.com
659    peaky = y;
660    if y <= npy-1 && y >= 1
661        f0 = log(result_conv(y,x));
662        f1 = real(log(result_conv(y-1,x)));
663        f2 = real(log(result_conv(y+1,x)));
664        peaky = peaky+ (f1-f2)/(2*f1-4*f0+2*f2);
665    else
666        F=-2; % warning flag for vector truncated by the limited search box
667    end
668    peakx=x;
669    if x <= npx-1 && x >= 1
670        f0 = log(result_conv(y,x));
671        f1 = real(log(result_conv(y,x-1)));
672        f2 = real(log(result_conv(y,x+1)));
673        peakx = peakx+ (f1-f2)/(2*f1-4*f0+2*f2);
674    else
675        F=-2; % warning flag for vector truncated by the limited search box
676    end
677    vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
678
679%------------------------------------------------------------------------
680% --- Find the maximum of the correlation function after interpolation
681function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
682%------------------------------------------------------------------------
683vector=[0 0]; %default
684F=-2;
685peaky=y;
686peakx=x;
687[npy,npx]=size(result_conv);
688if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1)
689    F=0;
690    for i=-1:1
691        for j=-1:1
692            %following 15 lines based on
693            %H. Nobach ï¿œ M. Honkanen (2005)
694            %Two-dimensional Gaussian regression for sub-pixel displacement
695            %estimation in particle image velocimetry or particle position
696            %estimation in particle tracking velocimetry
697            %Experiments in Fluids (2005) 38: 511ï¿œ515
698            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
699            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
700            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
701            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
702            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
703        end
704    end
705    c10=(1/6)*sum(sum(c10));
706    c01=(1/6)*sum(sum(c01));
707    c11=(1/4)*sum(sum(c11));
708    c20=(1/6)*sum(sum(c20));
709    c02=(1/6)*sum(sum(c02));
710    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
711    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
712    if abs(deltax)<1
713        peakx=x+deltax;
714    end
715    if abs(deltay)<1
716        peaky=y+deltay;
717    end
718end
719vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
720
721%'RUN_FIX': function for fixing velocity fields:
722%-----------------------------------------------
723% RUN_FIX(filename,field,flagindex,thresh_vecC,thresh_vel,iter,flag_mask,maskname,fileref,fieldref)
724%
725%filename: name of the netcdf file (used as input and output)
726%field: structure specifying the names of the fields to fix (depending on civ1 or civ2)
727    %.vel_type='civ1' or 'civ2';
728    %.nb=name of the dimension common to the field to fix ('nb_vectors' for civ1);
729    %.fixflag=name of fix flag variable ('vec_FixFlag' for civ1)
730%flagindex: flag specifying which values of vec_f are removed:
731        % if flagindex(1)=1: vec_f=-2 vectors are removed
732        % if flagindex(2)=1: vec_f=3 vectors are removed
733        % if flagindex(3)=1: vec_f=2 vectors are removed (if iter=1) or vec_f=4 vectors are removed (if iter=2)
734%iter=1 for civ1 fields and iter=2 for civ2 fields
735%thresh_vecC: threshold in the image correlation vec_C
736%flag_mask: =1 mask used to remove vectors (0 else)
737%maskname: name of the mask image file for fix
738%thresh_vel: threshold on velocity, or on the difference with the reference file fileref if exists
739%inf_sup=1: remove values smaller than threshold thresh_vel, =2, larger than threshold
740%fileref: .nc file name for a reference velocity (='': refrence 0 used)
741%fieldref: 'civ1','filter1'...feld used in fileref
742
743function FF=fix(Param,F,C,U,V,X,Y)
744FF=zeros(size(F));%default
745
746%criterium on warn flags
747FlagName={'CheckFmin2','CheckF2','CheckF3','CheckF4'};
748FlagVal=[-2 2 3 4];
749for iflag=1:numel(FlagName)
750    if isfield(Param,FlagName{iflag}) && Param.(FlagName{iflag})
751        FF=(FF==1| F==FlagVal(iflag));
752    end
753end
754%criterium on correlation values
755if isfield (Param,'MinCorr')
756    FF=FF==1 | C<Param.MinCorr;
757end
758if (isfield(Param,'MinVel')&&~isempty(Param.MinVel))||(isfield (Param,'MaxVel')&&~isempty(Param.MaxVel))
759    Umod= U.*U+V.*V;
760    if isfield (Param,'MinVel')&&~isempty(Param.MinVel)
761        FF=FF==1 | Umod<(Param.MinVel*Param.MinVel);
762    end
763    if isfield (Param,'MaxVel')&&~isempty(Param.MaxVel)
764        FF=FF==1 | Umod>(Param.MaxVel*Param.MaxVel);
765    end
766end
767
768
769
770
771
772
773
Note: See TracBrowser for help on using the repository browser.