source: trunk/src/civ_matlab.m @ 526

Last change on this file since 526 was 497, checked in by sommeria, 12 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.