source: trunk/src/civ_matlab.m @ 497

Last change on this file since 497 was 497, checked in by sommeria, 9 years ago

cleaning and small bug repair.
pb of histogram for filter data solved
display of uicontrol by right mouse selection improved

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