source: trunk/src/civ_matlab.m @ 837

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