source: trunk/src/civ_matlab.m @ 356

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

civ updated with new functions for opening files, consistently with uvmat
Bugs to be expected (use previous version then)

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