source: trunk/src/civ_matlab.m @ 494

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

various bugs corrected after testing in Windows OS. Introduction
of filter tps

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