YaDICs  V04.14.a
Yet another Digital Image Correlation software: platform dedicated to 2/3D Fluid and Solid kinematics field measurements.
 All Classes Files Functions Variables Pages
Ccorrelation_opticalFlow.h
Go to the documentation of this file.
1 /**********************************************************************
2  * Copyright (C) 2012, The YaDICs Project Developers.
3  * See the COPYRIGHT file at the top-level directory of this distribution ./COPYRIGHT.
4  * See ./COPYING file for copying and redistribution conditions.
5  *
6  * This file is part of YaDICs.
7  * YaDICs is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * (at your option) any later version.
11  *
12  * YaDICs is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with YaDICs. If not, see <http://www.gnu.org/licenses/>.
19  *
20  * Information about how to use the software are provided at http://yadic.univ-lille1.fr/
21  **********************************************************************/
22 
23 
24 #ifndef CCORRELATION_OPTICALFLOW
25 #define CCORRELATION_OPTICALFLOW
26 
27 
37 #include "Ccorrelation.h"
38 
52 template<typename T,typename Timg>
54 {
55 
56  public:
57 
58  std::string m_shape_type, m_interpLoop, m_interpFinal;
59  double m_epsilon;
60  std::vector<T> m_ref;
61 
63  {// constructor
64 
65  this->class_name = "Correlation : Optical Flow";
66 
67  }
68 
69  virtual void init(CParameterNetCDF &fp)
70  {// init from parameter file
71 
76  std::string attribute_name("shape");
77  int error = fp.loadAttribute(attribute_name,m_shape_type);
78  if (this->m_verbose){printf("\tShape function type : %s\n", m_shape_type.c_str());}
79 
80  attribute_name = "convergence_speed";
81  error = fp.loadAttribute(attribute_name,m_epsilon);
82  if (this->m_verbose){printf("\tLimit speed of convergence : %f\n", m_epsilon);}
83 
84  attribute_name = "interpolation_loop";
85  m_interpLoop = "linear";
86  error = fp.loadAttribute(attribute_name,m_interpLoop);
87  if (this->m_verbose){printf("\tLoop interpolator : %f\n", m_interpLoop.c_str());}
88 
89  attribute_name = "interpolation_final";
90  m_interpFinal = "linear";
91  error = fp.loadAttribute(attribute_name,m_interpFinal);
92  if (this->m_verbose){printf("\tFinal iterpolator : %f\n", m_interpFinal.c_str());}
93 
94 
95  }
96 
97  virtual void exec(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, Cfield<T,Timg> &oField)
98  {// exec optical flow correlation
99 
100  //-----------------------------------------------------------------------
101  //monoprocessing asignments
102  this->_3D = oMesh._3D;
103 
104  int iter(0); //number of iterations
105  int converged(0); //convergence parameter
106  T previous_res(1e10); //residue at previous step
107 
108  CImgList<T> K,Ke;//global and local tensors for linear system K[0].U=K[1] with K[0] = sum_(V.Vt) and K[1] = sum_((Img1(X)-Img2(X-U)).V)
109  CImg<int> p(this->m_threadNB,1,1,1,0);
110  CImgList<T> solution(oField.m_mode);//intermediate solution vector
111  CImg<T> current_res(this->m_threadNB,1,1,1,0);//global mean residue
112  assignG(oMesh, K,Ke);//global assignments
113  //-----------------------------------------------------------------------
114 
115  #pragma omp parallel
116  {//beginning of parallel section
117 
118  //-----------------------------------------------------------------------
119  //multithread initializations
120  //for elements
121  CImg<Timg> list; //elementary pixel list
122 
123  //for pixels
124  CImgList<T> Ve; //shape function to image gradient product at pix
125  CImgList<T> N(oMesh.num_var(), oMesh.num_mod(),1,1,1, 0);//Shape function at pix
126  T residue(0), deformed(0); //residue and corrected GL at pix
127  std::vector<T> ref, offset, lCoor, gCoor, grad, U;
128  ref.assign(3,0);offset.assign(3,0);lCoor.assign(3,0); gCoor.assign(3,0); U.assign(3,0); grad.assign(oMesh.num_var(),0);
129 
130  assignL(oMesh, Ve);
131  //-----------------------------------------------------------------------
132 
133  //calcul init value of the metric
134  metric(oImage,oMesh,solution,oField,N,Ve,Ke,ref,offset,lCoor,gCoor,U,grad,residue,deformed,list, p,current_res,K);
135 
136 
137  do
138  {//while the convergency velocity is higher than a certain level
139 
140  #pragma omp single
141  {
142  solve(K,oMesh,solution);
143  if (this->m_medianFilter){this->regularize(oMesh, oImage, solution);}
144  }
145 
146  metric(oImage,oMesh,solution,oField,N,Ve,Ke,ref,offset,lCoor,gCoor,U,grad,residue,deformed,list, p,current_res,K);
147 
148  #pragma omp single
149  {
150  test(solution, m_epsilon, oField, previous_res, iter, converged);
151  }
152 
153  }
154  while(!converged);
155 
156  if (oImage.m_storeActiv)
157  {export_res(oImage,oMesh,oField, N,ref,offset,lCoor,gCoor,U,residue,deformed,list);}
158 
159  }// pragma omp parallel
160 
161  }
162 
163  virtual void assignG(const Cmesh<T,Timg> &oMesh, CImgList<T> &K, CImgList<T> &Ke)
164  {//assign approriate size for Ke and Fe tensors
165 
166  K.assign(2);
167  Ke.assign(2);
168 
169  }
170 
171  virtual void assignL(const Cmesh<T,Timg> &oMesh, CImgList<T> &Ve) = 0;
172 
173  virtual void reset(const int &thread, CImgList<T> &Ke)
174  {//reset tensors
175 
176  cimglist_for(Ke, var)
177  {
178  cimg_forXYZ(Ke[var], x,y,z)
179  {
180  Ke[var](x,y,z,thread) = 0;
181  }
182  }
183 
184  }
185 
186  virtual void evalGN(const std::vector<T> &grad,const CImgList<T> &N, CImgList<T> &Ve) = 0;
187 
188  virtual void e_ref(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, const int &elem, std::vector<T> &ref) = 0;
189 
190  virtual void evalU(const CImgList<T> &solution, const Cmesh<T,Timg> &oMesh, const int &elem, const CImgList<T> &N, std::vector<T> &U) = 0;
191 
192  virtual void evalShape(const CImg<Timg> &list, const Cimage<Timg> &oImage, const Cfield<T,Timg> &oField, const std::vector<T> &ref, const std::vector<T> &lCoor, const std::vector<T> &gCoor, CImgList<T> &N) = 0;
193 
194  virtual void gradient(const Cimage<Timg> &oImage, const std::vector<T> &x, std::vector<T> &grad)
195  {// image gradient calcul at "x, y, z" image location
196 
201  grad[0] = (T)( oImage.m_curImg[0](x[0]+1,x[1],x[2], 0) - oImage.m_curImg[0](x[0]-1,x[1],x[2], 0) )/(2.*oImage.m_pix[0]);
202  grad[1] = (T)( oImage.m_curImg[0](x[0],x[1]+1,x[2], 0) - oImage.m_curImg[0](x[0],x[1]-1,x[2], 0) )/(2.*oImage.m_pix[1]);
203 
204  if (this->_3D)
205  {grad[2] = (T)( oImage.m_curImg[0](x[0],x[1],x[2]+1, 0) - oImage.m_curImg[0](x[0],x[1],x[2]-1, 0) )/(2.*oImage.m_pix[2]);}
206 
207  }
208 
209  void evalRes_loop(const Cimage<Timg> &oImage,const std::vector<T> &U, const std::vector<T> &x, T &residue, T &deformed)
210  {//eval Corrected deformed image and residue at pix
211 
212  //corrected deformed image
214  if (m_interpLoop == "cubic")
215  {
216  deformed = oImage.m_curImg[1].cubic_atXYZ(x[0]+U[0]/oImage.m_pix[0],x[1]+U[1]/oImage.m_pix[1],x[2]+U[2]/oImage.m_pix[2], 0);
217  }
218  else
219  {
220  deformed = oImage.m_curImg[1].linear_atXYZC(x[0]+U[0]/oImage.m_pix[0],x[1]+U[1]/oImage.m_pix[1],x[2]+U[2]/oImage.m_pix[2], 0,0);
221  }
222 
223  //residue between ref and corrected deformed image
224  residue = 0;
225  if (deformed != 0)
226  {residue = (T)oImage.m_curImg(0, x[0],x[1],x[2])- deformed;}
227 
228  }
229 
230  void evalRes_final(const Cimage<Timg> &oImage,const std::vector<T> &U, const std::vector<T> &x, T &residue, T &deformed)
231  {//eval Corrected deformed image and residue at pix
232 
234  if (m_interpFinal == "cubic")
235  {
236  deformed = oImage.m_curImg[1].cubic_atXYZ(x[0]+U[0]/oImage.m_pix[0],x[1]+U[1]/oImage.m_pix[1],x[2]+U[2]/oImage.m_pix[2],0);
237  }
238  else
239  {
240  deformed = oImage.m_curImg[1].linear_atXYZC(x[0]+U[0]/oImage.m_pix[0],x[1]+U[1]/oImage.m_pix[1],x[2]+U[2]/oImage.m_pix[2], 0,0);
241  }
242 
243  //residue between ref and corrected deformed image
244  residue = 0;
245  if (deformed != 0)
246  {residue = (T)oImage.m_curImg(0, x[0],x[1],x[2])- deformed;}
247 
248  }
249 
250 
251  virtual void tensors(const int &thread, const CImgList<T> &Ve, const T &residue, CImgList<T> &Ke)
252  {//eval stiffness matrix & force vector at element
253 
254  cimg_forZ(Ke[0], dim)
255  {//in the case of OFI dim = 0
256  cimg_forX(Ke[0],modi)
257  {//for every modes or nodes
258  cimg_forY(Ke[0],modj)
259  {
260  Ke[0](modi,modj,dim,thread) += Ve[dim][modi] * Ve[dim][modj];//Ke
261  }
262  Ke[1](0, modi,dim,thread) += Ve[dim][modi] * residue;//Fe
263  }
264  }
265 
266  }
267 
268  virtual void assembly(const int &thread, const Cmesh<T,Timg> &oMesh, const int &elem, const CImgList<T> &Ke, CImgList<T> &K) = 0;
269 
270  virtual void metric(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, CImgList<T> &solution, Cfield<T,Timg> &oField, CImgList<T> &N, CImgList<T> &Ve, CImgList<T> &Ke, std::vector<T> &ref, std::vector<T> &offset, std::vector<T> &lCoor, std::vector<T> &gCoor, std::vector<T> &U, std::vector<T> &grad, T &residue, T &deformed, CImg<Timg> &list, CImg<int> &p, CImg<T> &current_res, CImgList<T> &K)
271  {//calcul K
272 
273  #pragma omp single
274  {//init
275  current_res.fill(0);
276  cimglist_for(K,var){K[var].fill(0);}
277  }
278 
279  int thread(0),enter(1);
280 
281  #pragma omp for schedule(runtime)
282  cimglist_for(oMesh.m_connect, elem)
283  {//for every elements
284 
285  if (oMesh.cropList(oImage, elem, list, offset))
286  {//if the list is not entirely masked
287 
288  //check the thread number
289  thread = omp_get_thread_num();
290 
291  //look for the reference system of the element
292  e_ref(oImage, oMesh, elem, ref);
293 
294  //reset tensors
295  reset(thread,Ke);
296 
297  for (int pix=0; pix<list.size(); pix++)
298  {//eval Ke over pix
299 
300  //define gloMesh.obal coordonates : find back "x,y,z" triplet
301  list.contains(list[pix],lCoor[0],lCoor[1],lCoor[2]);
302  for (int i=0; i<gCoor.size(); i++){gCoor[i]=lCoor[i]+offset[i];}
303 
304  enter = 1;
305  //if pix is unmasked
306  if (!oImage.m_mask.is_empty()){enter = oImage.m_curMask(gCoor[0], gCoor[1], gCoor[2]);}
307 
308  if (enter)
309  {//if pixels are unmasked and in bounds
310 
311  // Evaluates grad[dim] at (x,y,z) in physical unit
312  gradient(oImage, gCoor, grad);
313 
314  // Evaluates N[dim][mode] at (x,y,z) in physical unit
315  evalShape(list, oImage, oField, ref, lCoor, gCoor, N);
316 
317  // Evaluates Ve[dim][mode] at (x,y,z) in physical unit
318  evalGN(grad,N, Ve);
319 
320  // Evaluates U[dim] in physical unit
321  evalU(solution, oMesh, elem, N, U);
322 
323  // Evaluates residue between ref and corrected deformed image at (x,y,z)
324  evalRes_loop(oImage, U, gCoor, residue, deformed);
325 
326  //sum of square residues over pixels and elements
327  current_res[thread] += std::pow(residue,2);
328 
329  //evaluate Ke
330  tensors(thread, Ve, residue, Ke);
331 
332  }//endif unmasked or out-of-bound
333 
334  }//for implied barrier
335 
336  assembly(thread, oMesh, elem, Ke, K);
337 
338  }
339 
340  }
341 
342  #pragma omp single
343  {//sum of square residues over threads and normalization and store residue of previous increment
344  oField.m_meanRes.push_back(std::sqrt(current_res.sum())*100/oImage.m_res);
345  }
346 
347  }
348 
349  virtual void solve(CImgList<T> &K, const Cmesh<T,Timg> &oMesh, CImgList<T> &solution)
350  {// K.U = F
351 
352  solution = K[1].get_solve(K[0]);
353 
354  }
355 
356  virtual void test(const CImgList<T> &solution, const double &epsilon, Cfield<T,Timg> &oField, T &previous_res, int &iter, int &converged)
357  {//convergency test
358 
359  if (std::abs(oField.m_meanRes.back() - previous_res) >= epsilon and oField.m_meanRes.back()<=previous_res)
360  {
361 
362  iter++;
363  previous_res = oField.m_meanRes.back();
364  oField.m_mode = solution;
365 
366  if(this->m_verbose){std::cout<<"(iter "<<iter<<") : "<<oField.m_meanRes.back()<<" \%\n\n";}
367 
368  }
369  else
370  {
371 
372  converged = 1;
373  if (oField.m_meanRes.back()>previous_res){printf("Warning : last step has diverged\n");}
374  else{printf("Last step has converged\n");}
375  oField.m_meanRes.back() = previous_res;
376  std::cout<<this->class_name.c_str()<<" -> ultimate residue ("<<iter<<" iter) : "<<previous_res<<" \%\n";
377 
378  }
379 
380  }
381 
382  virtual void export_res(const Cimage<Timg> &oImage, const Cmesh<T,Timg> &oMesh, Cfield<T,Timg> &oField, CImgList<T> &N, std::vector<T> &ref, std::vector<T> &offset, std::vector<T> &lCoor, std::vector<T> &gCoor, std::vector<T> &U, T &residue, T &deformed, CImg<Timg> &list)
383  {//actualize residue
384 
385  int thread(0), enter(1);
386 
387  #pragma omp for schedule(runtime)
388  cimglist_for(oMesh.m_connect, elem)
389  {//for every elements
390 
391  if (oMesh.cropList(oImage, elem, list, offset))
392  {//if the list is not entirely masked
393 
394  //check the thread number
395  thread = omp_get_thread_num();
396 
397  //look for the reference system of the element
398  e_ref(oImage, oMesh, elem, ref);
399 
400  for (int pix=0; pix<list.size(); pix++)
401  {//eval Ke over pix
402 
403  //define global coordonates : find back "x,y,z" triplet
404  list.contains(list[pix],lCoor[0],lCoor[1],lCoor[2]);
405  for (int i=0; i<gCoor.size(); i++){gCoor[i]=lCoor[i]+offset[i];}
406 
407  enter = 1;
408  //if pix is unmasked
409  if (!oImage.m_mask.is_empty()){enter = oImage.m_curMask(gCoor[0], gCoor[1], gCoor[2]);}
410 
411  if (enter)
412  {//if pixels are unmasked and in bounds
413 
414  // Evaluates N[dim][mode] at (x,y,z) in physical unit
415  evalShape(list, oImage, oField, ref, lCoor, gCoor, N);
416 
417  // Evaluates U[dim] in physical unit
418  evalU(oField.m_mode, oMesh, elem, N, U);
419 
420  // Evaluates residue between ref and corrected deformed image at (x,y,z)
421  evalRes_final(oImage, U, gCoor, residue, deformed);
422 
423  //export data
424  store(oField,oImage,gCoor,U,deformed,residue);
425 
426  }
427 
428  }//for
429 
430  }//if croplist
431 
432  }//for
433 
434  }
435 
436  virtual void store(Cfield<T,Timg> &oField, const Cimage<Timg> &oImage, const std::vector<T> gCoor, const std::vector<T> &U,const T &deformed,const T &residue)
437  {//store output data
438 
439  if (oImage.m_storeField)
440  {//store field a iter - 1 to avoid storing the last diverging increment
441  cimglist_for(oField.m_field, dim)
442  {
443  #pragma omp critical
444  {
445  oField.m_field(dim, gCoor[0],gCoor[1],gCoor[2]) = U[dim];
446  }
447  }
448  }
449 
450  if (oImage.m_storeDef)
451  {//store imDef a iter - 1
452 // #pragma omp critical
453 // {
454  oField.m_def(gCoor[0],gCoor[1],gCoor[2]) = deformed;
455 // }
456  }
457 
458  if (oImage.m_storeRes)
459  {//store imRes a iter - 1
460 // #pragma omp critical
461 // {
462  oField.m_res(gCoor[0],gCoor[1],gCoor[2]) = residue;
463 // }
464  }
465 
466  }
467 
468 
469 };
470 
471 #endif
virtual void regularize(const Cmesh< T, Timg > &oMesh, const Cimage< Timg > &oImage, CImgList< T > &solution)
Definition: Ccorrelation.h:96
void evalRes_final(const Cimage< Timg > &oImage, const std::vector< T > &U, const std::vector< T > &x, T &residue, T &deformed)
Definition: Ccorrelation_opticalFlow.h:230
void evalRes_loop(const Cimage< Timg > &oImage, const std::vector< T > &U, const std::vector< T > &x, T &residue, T &deformed)
Definition: Ccorrelation_opticalFlow.h:209
this is the mother class of every correlation classes. m_correl_type refers to the the type of correl...
Definition: Ccorrelation.h:65
This is the mother class of Optical Flow strategies: It implements different steps required for the i...
Definition: Ccorrelation_opticalFlow.h:53
Definition: Cmesh.h:61
Definition: Cfield.h:51
Definition: Cimage.h:57
virtual void gradient(const Cimage< Timg > &oImage, const std::vector< T > &x, std::vector< T > &grad)
Definition: Ccorrelation_opticalFlow.h:194
virtual void init(CParameterNetCDF &fp)
Definition: Ccorrelation_opticalFlow.h:69