Visual Servoing Platform  version 3.2.1 under development (2019-05-26)
vpTemplateTrackerMIESM.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2019 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See http://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * Example of template tracking.
33  *
34  * Authors:
35  * Amaury Dame
36  * Aurelien Yol
37  * Fabien Spindler
38  *
39  *****************************************************************************/
40 
41 #include <visp3/tt_mi/vpTemplateTrackerMIESM.h>
42 
43 #ifdef VISP_HAVE_OPENMP
44 #include <omp.h>
45 #endif
46 
48  : vpTemplateTrackerMI(_warp), minimizationMethod(USE_NEWTON), CompoInitialised(false), HDirect(), HInverse(),
49  HdesireDirect(), HdesireInverse(), GDirect(), GInverse()
50 {
51  useCompositionnal = false;
52  useInverse = false;
53  if (!Warp->isESMcompatible())
54  std::cerr << "The selected warp function is not appropriate for the ESM "
55  "algorithm..."
56  << std::endl;
57 }
58 
60 {
62  std::cout << "Initialise Hessian at Desired position..." << std::endl;
63 
64  dW = 0;
65 
66  // double erreur=0;
67  int Nbpoint = 0;
68 
69  double i2, j2;
70  // double Tij;
71  double /*IW,*/ dx, dy;
72  // int cr,ct;
73  // double er,et;
74  int i, j;
75 
76  Nbpoint = 0;
77  // erreur=0;
78 
79  if (blur)
81 
83 
84  Warp->computeCoeff(p);
85  for (unsigned int point = 0; point < templateSize; point++) {
86  i = ptTemplate[point].y;
87  j = ptTemplate[point].x;
88  X1[0] = j;
89  X1[1] = i;
90 
91  Warp->computeDenom(X1, p);
92  Warp->warpX(X1, X2, p);
93 
94  j2 = X2[0];
95  i2 = X2[1];
96 
97  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
98  Nbpoint++;
99  // if(blur)
100  // IW=BI.getValue(i2,j2);
101  // else
102  // IW=I.getValue(i2,j2);
103 
104  // ct=ptTemplateSupp[point].ct;
105  // et=ptTemplateSupp[point].et;
106  // cr=(int)((IW*(Nc-1))/255.);
107  // er=((double)IW*(Nc-1))/255.-cr;
108 
109  // if(ApproxHessian==vpTemplateTrackerMI::HESSIAN_NONSECOND)
110  // // cas à tester AY
111  // vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout,cr,
112  // er, ct, et, Nc, ptTemplate[point].dW, nbParam,
113  // bspline);
114  // else
115  // vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout,cr,
116  // er, ct, et, Nc, ptTemplate[point].dW, nbParam,
117  // bspline);
118 
119  // vpTemplateTrackerMIBSpline::computeProbabilities(PrtTout,cr,
120  // er, ct, et, Nc, ptTemplate[point].dW,
121  // nbParam,bspline,ApproxHessian,false);
122  }
123  }
124 
125  double MI;
126  computeProba(Nbpoint);
127  computeMI(MI);
129  // std::cout<<"HdesireInverse"<<std::endl<<HdesireInverse<<std::endl;
130 
132  // DIRECT COMPO
133 
141  }
142 
143  Nbpoint = 0;
144  // erreur=0;
146 
147  Warp->computeCoeff(p);
148  for (unsigned int point = 0; point < templateSize; point++) {
149  i = ptTemplate[point].y;
150  j = ptTemplate[point].x;
151  X1[0] = j;
152  X1[1] = i;
153 
154  Warp->computeDenom(X1, p);
155  Warp->warpX(X1, X2, p);
156 
157  j2 = X2[0];
158  i2 = X2[1];
159 
160  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight()) && (j2 < I.getWidth())) {
161  Nbpoint++;
162  // Tij=ptTemplate[point].val;
163  // if(!blur)
164  // IW=I.getValue(i2,j2);
165  // else
166  // IW=BI.getValue(i2,j2);
167 
168  dx = 1. * dIx.getValue(i2, j2) * (Nc - 1) / 255.;
169  dy = 1. * dIy.getValue(i2, j2) * (Nc - 1) / 255.;
170 
171  // cr=ptTemplateSupp[point].ct;
172  // er=ptTemplateSupp[point].et;
173  // ct=(int)((IW*(Nc-1))/255.);
174  // et=((double)IW*(Nc-1))/255.-ct;
175 
176  Warp->dWarpCompo(X1, X2, p, ptTemplateCompo[point].dW, dW);
177 
178  double *tptemp = new double[nbParam];
179  for (unsigned int it = 0; it < nbParam; it++)
180  tptemp[it] = dW[0][it] * dx + dW[1][it] * dy;
181 
182  // vpTemplateTrackerMIBSpline::computeProbabilities(PrtTout,cr,
183  // er, ct, et, Nc, ptTemplate[point].dW,
184  // nbParam,bspline,ApproxHessian,
185  // hessianComputation==
186  // vpTemplateTrackerMI::USE_HESSIEN_DESIRE);
187  // calcul de l'erreur
188  // erreur+=(Tij-IW)*(Tij-IW);
189 
190  // if(ApproxHessian==vpTemplateTrackerMI::HESSIAN_NONSECOND) //
191  // cas à tester AY
192  // vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout,cr,
193  // er, ct, et, Nc, tptemp, nbParam, bspline);
194  // else
195  // vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout,cr,
196  // er, ct, et, Nc, tptemp, nbParam, bspline);
197 
198  // vpTemplateTrackerMIBSpline::computeProbabilities(PrtTout,cr, er,
199  // ct, et, Nc, tptemp, nbParam,bspline,ApproxHessian,false);
200 
201  delete[] tptemp;
202  }
203  }
204 
205  computeProba(Nbpoint);
206  computeMI(MI);
208  // std::cout<<"HdesireDirect"<<std::endl<<HdesireDirect<<std::endl;
209 
210  lambda = lambdaDep;
211 
213  // Hdesire=HdesireDirect;
214  // Hdesire=HdesireInverse;
215 
216  // std::cout<<"HdesireDirect "<<HdesireDirect<<std::endl;
217  // std::cout<<"HdesireInverse "<<HdesireInverse<<std::endl;
220 
221  // std::cout<<"\tEnd initialisation..."<<std::endl;
222 }
223 
225 {
226  // std::cout<<"Initialise precomputed value of ESM Tracker"<<std::endl;
227 
234 
235  ptTemplateSupp = new vpTemplateTrackerPointSuppMIInv[templateSize];
237  for (unsigned int point = 0; point < templateSize; point++) {
238  int i = ptTemplate[point].y;
239  int j = ptTemplate[point].x;
240  X1[0] = j;
241  X1[1] = i;
242  Warp->computeDenom(X1, p);
243 
244  ptTemplateCompo[point].dW = new double[2 * nbParam];
245  Warp->getdWdp0(i, j, ptTemplateCompo[point].dW);
246 
247  ptTemplate[point].dW = new double[nbParam];
248  double dx = ptTemplate[point].dx * (Nc - 1) / 255.;
249  double dy = ptTemplate[point].dy * (Nc - 1) / 255.;
250  Warp->getdW0(i, j, dy, dx, ptTemplate[point].dW);
251 
252  double Tij = ptTemplate[point].val;
253  int ct = (int)((Tij * (Nc - 1)) / 255.);
254  double et = (Tij * (Nc - 1)) / 255. - ct;
255  ptTemplateSupp[point].et = et;
256  ptTemplateSupp[point].ct = ct;
257  ptTemplateSupp[point].Bt = new double[4];
258  ptTemplateSupp[point].dBt = new double[4];
259  for (char it = -1; it <= 2; it++) {
260  ptTemplateSupp[point].Bt[it + 1] = vpTemplateTrackerBSpline::Bspline4(-it + et);
261  ptTemplateSupp[point].dBt[it + 1] = vpTemplateTrackerMIBSpline::dBspline4(-it + et);
262  }
263  }
264  CompoInitialised = true;
265 }
266 
268 {
269  if (!CompoInitialised)
270  std::cout << "Compositionnal tracking no initialised\nUse "
271  "initCompInverse(vpImage<unsigned char> &I) function"
272  << std::endl;
273  dW = 0;
274 
275  if (blur)
279  /* if(ApproxHessian!=HESSIAN_NONSECOND && ApproxHessian!=HESSIAN_0 &&
280  ApproxHessian!=HESSIAN_NEW && ApproxHessian!=HESSIAN_YOUCEF)
281  {
282  getGradX(dIx, d2Ix,fgdG,taillef);
283  getGradY(dIx, d2Ixy,fgdG,taillef);
284  getGradY(dIy, d2Iy,fgdG,taillef);
285  }*/
286 
287  // double erreur=0;
288  int point;
289 
290  MI_preEstimation = -getCost(I, p);
291 
292  lambda = lambdaDep;
293  // double MIprec=-1000;
294 
295  double i2, j2;
296  // double Tij;
297  // double IW;
298  // int cr,ct;
299  // double er,et;
300 
301  vpColVector dpinv(nbParam);
302 
303  double alpha = 2.;
304 
305  int i, j;
306  unsigned int iteration = 0;
307  do {
308  int Nbpoint = 0;
309  // MIprec=MI;
310  double MI = 0;
311  // erreur=0;
312 
314 
316  // Inverse
317  Warp->computeCoeff(p);
318  for (point = 0; point < (int)templateSize; point++) {
319  i = ptTemplate[point].y;
320  j = ptTemplate[point].x;
321  X1[0] = j;
322  X1[1] = i;
323 
324  Warp->computeDenom(X1, p);
325  Warp->warpX(X1, X2, p);
326 
327  j2 = X2[0];
328  i2 = X2[1];
329 
330  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
331  Nbpoint++;
332  // Tij=ptTemplate[point].val;
333  // if(!blur)
334  // IW=I.getValue(i2,j2);
335  // else
336  // IW=BI.getValue(i2,j2);
337 
338  // ct=ptTemplateSupp[point].ct;
339  // et=ptTemplateSupp[point].et;
340  // cr=(int)((IW*(Nc-1))/255.);
341  // er=((double)IW*(Nc-1))/255.-cr;
342 
343  // if(ApproxHessian==vpTemplateTrackerMI::HESSIAN_NONSECOND||hessianComputation==vpTemplateTrackerMI::USE_HESSIEN_DESIRE)
344  // // cas à tester AY
345  // vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout,
346  // cr, er, ct, et, Nc, ptTemplate[point].dW,
347  // nbParam, bspline);
348  // else
349  // vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout,
350  // cr, er, ct, et, Nc, ptTemplate[point].dW,
351  // nbParam, bspline);
352  }
353  }
354 
355  if (Nbpoint == 0) {
356  // std::cout<<"plus de point dans template suivi"<<std::endl;
357  diverge = true;
358  MI = 0;
359  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
360  } else {
361  computeProba(Nbpoint);
362  computeMI(MI);
365  computeGradient();
366  GInverse = G;
367 
369  // DIRECT
370 
371  Nbpoint = 0;
372  // MIprec=MI;
373  MI = 0;
374  // erreur=0;
375 
377 
378  Warp->computeCoeff(p);
379 #ifdef VISP_HAVE_OPENMP
380  int nthreads = omp_get_num_procs();
381  // std::cout << "file: " __FILE__ << " line: " << __LINE__ << "
382  // function: " << __FUNCTION__ << " nthread: " << nthreads << std::endl;
383  omp_set_num_threads(nthreads);
384 #pragma omp parallel for private(point, i, j, i2, j2) default(shared)
385 #endif
386  for (point = 0; point < (int)templateSize; point++) {
387  i = ptTemplate[point].y;
388  j = ptTemplate[point].x;
389  X1[0] = j;
390  X1[1] = i;
391  Warp->computeDenom(X1, p);
392  Warp->warpX(i, j, i2, j2, p);
393  X2[0] = j2;
394  X2[1] = i2;
395 
396  // Warp->computeDenom(X1,p);
397  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
398  Nbpoint++;
399  // Tij=ptTemplate[point].val;
400  // Tij=Iterateurvecteur->val;
401  // if(!blur)
402  // IW=I.getValue(i2,j2);
403  // else
404  // IW=BI.getValue(i2,j2);
405 
406  double dx = 1. * dIx.getValue(i2, j2) * (Nc - 1) / 255.;
407  double dy = 1. * dIy.getValue(i2, j2) * (Nc - 1) / 255.;
408 
409  // ct=(int)((IW*(Nc-1))/255.);
410  // et=((double)IW*(Nc-1))/255.-ct;
411  // cr=ptTemplateSupp[point].ct;
412  // er=ptTemplateSupp[point].et;
413 
414  Warp->dWarpCompo(X1, X2, p, ptTemplateCompo[point].dW, dW);
415 
416  double *tptemp = new double[nbParam];
417  for (unsigned int it = 0; it < nbParam; it++)
418  tptemp[it] = dW[0][it] * dx + dW[1][it] * dy;
419 
420  // calcul de l'erreur
421  // erreur+=(Tij-IW)*(Tij-IW);
422  // if(bspline==3)
423  // vpTemplateTrackerMIBSpline::PutTotPVBspline3NoSecond(PrtTout,
424  // cr, er, ct, et, Nc, tptemp, nbParam);
425  // else
426  // vpTemplateTrackerMIBSpline::PutTotPVBspline4NoSecond(PrtTout,
427  // cr, er, ct, et, Nc, tptemp, nbParam);
428  // vpTemplateTrackerMIBSpline::computeProbabilities(PrtTout,cr,
429  // er, ct, et, Nc, tptemp, nbParam,bspline,ApproxHessian,
430  // hessianComputation==vpHessienType::USE_HESSIEN_DESIRE);
431 
432  delete[] tptemp;
433  }
434  }
435 
436  computeProba(Nbpoint);
437  computeMI(MI);
440  computeGradient();
441  GDirect = G;
442 
444  H = HDirect + HInverse;
446  }
447  G = GDirect - GInverse;
448  // G=GDirect;
449  // G=-GInverse;
450 
451  try {
452  if (minimizationMethod == vpTemplateTrackerMIESM::USE_GRADIENT)
453  dp = -gain * 0.3 * G;
454  else {
455  switch (hessianComputation) {
457  dp = gain * HLMdesireInverse * G;
458  break;
460  if (HLM.cond() > HLMdesire.cond())
461  dp = gain * HLMdesireInverse * G;
462  else
463  dp = gain * 0.3 * HLM.inverseByLU() * G;
464  break;
465  default:
466  dp = gain * 0.3 * HLM.inverseByLU() * G;
467  break;
468  }
469  }
470  } catch (const vpException &e) {
471  // std::cerr<<"probleme inversion"<<std::endl;
472  throw(e);
473  }
474 
476  dp = -1. * dp;
477 
478  if (useBrent) {
479  alpha = 2.;
480  computeOptimalBrentGain(I, p, -MI, dp, alpha);
481  dp = alpha * dp;
482  }
483  p += dp;
484 
485  iteration++;
486  }
487  } while (/*(MI!=MIprec) &&*/ (iteration < iterationMax));
488 
489  MI_postEstimation = -getCost(I, p);
491  MI_postEstimation = -1;
492  }
493 
494  nbIteration = iteration;
495 }
virtual void dWarpCompo(const vpColVector &X1, const vpColVector &X2, const vpColVector &ParamM, const double *dwdp0, vpMatrix &dW)=0
void computeHessien(vpMatrix &H)
unsigned int getWidth() const
Definition: vpImage.h:237
static void getGradX(const vpImage< unsigned char > &I, vpImage< double > &dIx)
vpTemplateTrackerPoint * ptTemplate
static void getGradY(const vpImage< unsigned char > &I, vpImage< double > &dIy)
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
Definition: vpArray2D.h:305
vpTemplateTrackerMIESM()
Default constructor.
void computeOptimalBrentGain(const vpImage< unsigned char > &I, vpColVector &tp, double tMI, vpColVector &direction, double &alpha)
error that can be emited by ViSP classes.
Definition: vpException.h:71
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
static void getGradXGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIx, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
vpImage< double > d2Ix
vpImage< double > BI
vpHessienApproximationType ApproxHessian
unsigned int templateSize
unsigned int iterationMax
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1402
Error that can be emited by the vpTracker class and its derivates.
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp)
void initHessienDesired(const vpImage< unsigned char > &I)
vpImage< double > dIx
void computeMI(double &MI)
vpImage< double > dIy
vpImage< double > d2Iy
double cond(double svThreshold=1e-6) const
Definition: vpMatrix.cpp:5140
virtual bool isESMcompatible() const =0
virtual void getdW0(const int &i, const int &j, const double &dy, const double &dx, double *dIdW)=0
unsigned int nbIteration
void trackNoPyr(const vpImage< unsigned char > &I)
vpMinimizationTypeMIESM minimizationMethod
vpTemplateTrackerPointCompo * ptTemplateCompo
Implementation of column vector and the associated operations.
Definition: vpColVector.h:130
vpImage< double > d2Ixy
vpMatrix inverseByLU() const
unsigned int getHeight() const
Definition: vpImage.h:179
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, const bool convolve=false)
virtual void getdWdp0(const int &i, const int &j, double *dIdW)=0
vpTemplateTrackerWarp * Warp
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:5200
void computeProba(int &nbpoint)
vpHessienType hessianComputation
void resize(const unsigned int i, const bool flagNullify=true)
Definition: vpColVector.h:310