Visual Servoing Platform  version 3.2.0 under development (2019-01-22)
vpTemplateTrackerMIForwardAdditional.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/vpTemplateTrackerMIForwardAdditional.h>
42 
43 #ifdef VISP_HAVE_OPENMP
44 #include <omp.h>
45 #endif
46 
48  : vpTemplateTrackerMI(_warp), minimizationMethod(USE_NEWTON), evolRMS(0), x_pos(NULL), y_pos(NULL), threshold_RMS(0),
49  p_prec(), G_prec(), KQuasiNewton()
50 {
51  useCompositionnal = false;
52  threshold_RMS = 1e-20;
53 }
54 
56 {
57  // std::cout<<"Initialise Hessian at Desired position..."<<std::endl;
58 
59  dW = 0;
60 
61  int Nbpoint = 0;
62 
63  if (blur)
67 
68  double Tij;
69  double IW, dx, dy;
70  int cr, ct;
71  double er, et;
72 
73  Nbpoint = 0;
74 
76  Warp->computeCoeff(p);
77  for (unsigned int point = 0; point < templateSize; point++) {
78  int i = ptTemplate[point].y;
79  int j = ptTemplate[point].x;
80  X1[0] = j;
81  X1[1] = i;
82  X2[0] = j;
83  X2[1] = i;
84 
85  Warp->computeDenom(X1, p);
86  Warp->warpX(X1, X2, p);
87 
88  double j2 = X2[0];
89  double i2 = X2[1];
90 
91  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
92  Nbpoint++;
93  Tij = ptTemplate[point].val;
94  if (!blur)
95  IW = I.getValue(i2, j2);
96  else
97  IW = BI.getValue(i2, j2);
98 
99  dx = 1. * dIx.getValue(i2, j2) * (Nc - 1) / 255.;
100  dy = 1. * dIy.getValue(i2, j2) * (Nc - 1) / 255.;
101 
102  ct = (int)((IW * (Nc - 1)) / 255.);
103  cr = (int)((Tij * (Nc - 1)) / 255.);
104  et = (IW * (Nc - 1)) / 255. - ct;
105  er = ((double)Tij * (Nc - 1)) / 255. - cr;
106  // std::cout<<"test"<<std::endl;
107  Warp->dWarp(X1, X2, p, dW);
108 
109  double *tptemp = new double[nbParam];
110  for (unsigned int it = 0; it < nbParam; it++)
111  tptemp[it] = dW[0][it] * dx + dW[1][it] * dy;
112 
114  vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
116  vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
117 
118  delete[] tptemp;
119  }
120  }
121 
122  if (Nbpoint > 0) {
123  double MI;
124  computeProba(Nbpoint);
125  computeMI(MI);
127 
128  // double conditionnement=GetConditionnement(Hdesire);
129  // std::cout<<"conditionnement : "<<conditionnement<<std::endl;
131  try {
133  } catch (const vpException &e) {
134  // std::cerr<<"probleme inversion"<<std::endl;
135  throw(e);
136  }
137  // std::cout<<"Hdesire = "<<Hdesire<<std::endl;
138  // std::cout<<"\tEnd initialisation..."<<std::endl;
139  }
140 }
141 
143 {
144  dW = 0;
145 
146  // double erreur=0;
147  int Nbpoint = 0;
148  if (blur)
152 
153  double MI = 0, MIprec = -1000;
154 
155  MI_preEstimation = -getCost(I, p);
156 
157  double alpha = 2.;
158 
159  unsigned int iteration = 0;
160 
161  initPosEvalRMS(p);
162  do {
163  if (iteration % 5 == 0)
165  Nbpoint = 0;
166  MIprec = MI;
167  MI = 0;
168  // erreur=0;
169 
171 
172  Warp->computeCoeff(p);
173 #ifdef VISP_HAVE_OPENMP
174  int nthreads = omp_get_num_procs();
175  // std::cout << "file: " __FILE__ << " line: " << __LINE__ << " function:
176  // " << __FUNCTION__ << " nthread: " << nthreads << std::endl;
177  omp_set_num_threads(nthreads);
178 #pragma omp parallel for default(shared)
179 #endif
180  for (int point = 0; point < (int)templateSize; point++) {
181  int i = ptTemplate[point].y;
182  int j = ptTemplate[point].x;
183  X1[0] = j;
184  X1[1] = i;
185 
186  Warp->computeDenom(X1, p);
187  Warp->warpX(X1, X2, p);
188 
189  double j2 = X2[0];
190  double i2 = X2[1];
191 
192  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
193  Nbpoint++;
194  double Tij = ptTemplate[point].val;
195  double IW;
196  if (!blur)
197  IW = I.getValue(i2, j2);
198  else
199  IW = BI.getValue(i2, j2);
200 
201  double dx = 1. * dIx.getValue(i2, j2) * (Nc - 1) / 255.;
202  double dy = 1. * dIy.getValue(i2, j2) * (Nc - 1) / 255.;
203 
204  int ct = (int)((IW * (Nc - 1)) / 255.);
205  int cr = (int)((Tij * (Nc - 1)) / 255.);
206  double et = (IW * (Nc - 1)) / 255. - ct;
207  double er = ((double)Tij * (Nc - 1)) / 255. - cr;
208 
209  // calcul de l'erreur
210  // erreur+=(Tij-IW)*(Tij-IW);
211 
212  // Calcul de l'histogramme joint par interpolation bilinÈaire
213  // (Bspline ordre 1)
214  Warp->dWarp(X1, X2, p, dW);
215 
216  // double *tptemp=temp;
217  double *tptemp = new double[nbParam];
218  ;
219  for (unsigned int it = 0; it < nbParam; it++)
220  tptemp[it] = (dW[0][it] * dx + dW[1][it] * dy);
221  //*tptemp++ =dW[0][it]*dIWx+dW[1][it]*dIWy;
222  // std::cout<<cr<<" "<<ct<<" ; ";
224  vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
226  vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
227 
228  delete[] tptemp;
229  }
230  }
231 
232  if (Nbpoint == 0) {
233  // std::cout<<"plus de point dans template suivi"<<std::endl;
234  diverge = true;
235  MI = 0;
237  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
238  } else {
239  computeProba(Nbpoint);
240  computeMI(MI);
241  // std::cout<<iteration<<"\tMI= "<<MI<<std::endl;
242  computeHessien(H);
243  computeGradient();
244 
246  try {
247  switch (hessianComputation) {
249  dp = gain * HLMdesireInverse * G;
250  break;
252  if (HLM.cond() > HLMdesire.cond())
253  dp = gain * HLMdesireInverse * G;
254  else
255  dp = gain * 0.2 * HLM.inverseByLU() * G;
256  break;
257  default:
258  dp = gain * 0.2 * HLM.inverseByLU() * G;
259  break;
260  }
261  } catch (const vpException &e) {
262  // std::cerr<<"probleme inversion"<<std::endl;
264  throw(e);
265  }
266  }
267 
268  switch (minimizationMethod) {
270  vpColVector p_test_LMA(nbParam);
272  p_test_LMA = p - 100000.1 * dp;
273  else
274  p_test_LMA = p + 1. * dp;
275  MI = -getCost(I, p);
276  double MI_LMA = -getCost(I, p_test_LMA);
277  if (MI_LMA > MI) {
278  p = p_test_LMA;
279  lambda = (lambda / 10. < 1e-6) ? lambda / 10. : 1e-6;
280  } else {
281  lambda = (lambda * 10. < 1e6) ? 1e6 : lambda * 10.;
282  }
283  } break;
285  dp = -gain * 6.0 * G;
286  if (useBrent) {
287  alpha = 2.;
288  computeOptimalBrentGain(I, p, -MI, dp, alpha);
289  dp = alpha * dp;
290  }
291  p += 1. * dp;
292  break;
293  }
294 
296  if (iterationGlobale != 0) {
297  vpColVector s_quasi = p - p_prec;
298  vpColVector y_quasi = G - G_prec;
299  double s_scal_y = s_quasi.t() * y_quasi;
300  // if(s_scal_y!=0)//BFGS
301  // KQuasiNewton=KQuasiNewton-(s_quasi*y_quasi.t()*KQuasiNewton+KQuasiNewton*y_quasi*s_quasi.t())/s_scal_y+(1.+y_quasi.t()*(KQuasiNewton*y_quasi)/s_scal_y)*s_quasi*s_quasi.t()/s_scal_y;
302  // if(s_scal_y!=0)//DFP
303  if (std::fabs(s_scal_y) > std::numeric_limits<double>::epsilon())
304  KQuasiNewton = KQuasiNewton + 0.001 * (s_quasi * s_quasi.t() / s_scal_y -
305  KQuasiNewton * y_quasi * y_quasi.t() * KQuasiNewton /
306  (y_quasi.t() * KQuasiNewton * y_quasi));
307  }
308  dp = -KQuasiNewton * G;
309  p_prec = p;
310  G_prec = G;
311  p -= 1.01 * dp;
312  } break;
313 
314  default: {
316  dp = -0.1 * dp;
317  if (useBrent) {
318  alpha = 2.;
319  computeOptimalBrentGain(I, p, -MI, dp, alpha);
320  // std::cout<<alpha<<std::endl;
321  dp = alpha * dp;
322  }
323 
324  p += 1. * dp;
325  break;
326  }
327  }
328 
329  computeEvalRMS(p);
330  iteration++;
332 
333  } while ((std::fabs(MI - MIprec) > std::fabs(MI) * std::numeric_limits<double>::epsilon()) &&
334  (iteration < iterationMax) && (evolRMS > threshold_RMS));
335  // while( (MI!=MIprec) &&(iteration< iterationMax)&&(evolRMS>threshold_RMS)
336  // );
337  if (Nbpoint == 0) {
338  // std::cout<<"plus de point dans template suivi"<<std::endl;
340  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
341  }
342 
343  nbIteration = iteration;
344  MI_postEstimation = -getCost(I, p);
346  MI_postEstimation = -1;
347  }
349 }
350 
352 {
353  unsigned int nb_corners = zoneTracked->getNbTriangle() * 3;
354  x_pos = new double[nb_corners];
355  y_pos = new double[nb_corners];
356 
357  Warp->computeCoeff(pw);
358  vpTemplateTrackerTriangle triangle;
359 
360  for (unsigned int i = 0; i < zoneTracked->getNbTriangle(); i++) {
361  zoneTracked->getTriangle(i, triangle);
362  for (unsigned int j = 0; j < 3; j++) {
363  triangle.getCorner(j, X1[0], X1[1]);
364 
365  Warp->computeDenom(X1, pw);
366  Warp->warpX(X1, X2, pw);
367  x_pos[i * 3 + j] = X2[0];
368  y_pos[i * 3 + j] = X2[1];
369  }
370  }
371 }
372 
374 {
375  unsigned int nb_corners = zoneTracked->getNbTriangle() * 3;
376 
377  Warp->computeCoeff(pw);
378  evolRMS = 0;
379  vpTemplateTrackerTriangle triangle;
380 
381  for (unsigned int i = 0; i < zoneTracked->getNbTriangle(); i++) {
382  zoneTracked->getTriangle(i, triangle);
383  for (unsigned int j = 0; j < 3; j++) {
384  triangle.getCorner(j, X1[0], X1[1]);
385 
386  Warp->computeDenom(X1, pw);
387  Warp->warpX(X1, X2, pw);
388  evolRMS += (x_pos[i * 3 + j] - X2[0]) * (x_pos[i * 3 + j] - X2[0]) +
389  (y_pos[i * 3 + j] - X2[1]) * (y_pos[i * 3 + j] - X2[1]);
390  x_pos[i * 3 + j] = X2[0];
391  y_pos[i * 3 + j] = X2[1];
392  }
393  }
394  evolRMS = evolRMS / nb_corners;
395 }
396 
398 {
399  delete[] x_pos;
400  delete[] y_pos;
401 }
void computeHessien(vpMatrix &H)
void getTriangle(unsigned int i, vpTemplateTrackerTriangle &T) const
unsigned int getWidth() const
Definition: vpImage.h:239
void trackNoPyr(const vpImage< unsigned char > &I)
vpTemplateTrackerPoint * ptTemplate
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
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
vpColVector getCorner(unsigned int i) const
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 > BI
vpHessienApproximationType ApproxHessian
double cond() const
Definition: vpMatrix.cpp:5045
unsigned int templateSize
unsigned int iterationMax
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1420
Error that can be emited by the vpTracker class and its derivates.
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp)
unsigned int getNbTriangle() const
vpImage< double > dIx
void computeMI(double &MI)
vpRowVector t() const
vpImage< double > dIy
unsigned int iterationGlobale
unsigned int nbIteration
Implementation of column vector and the associated operations.
Definition: vpColVector.h:72
vpMatrix inverseByLU() const
vpTemplateTrackerZone * zoneTracked
unsigned int getHeight() const
Definition: vpImage.h:178
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, const bool convolve=false)
void initHessienDesired(const vpImage< unsigned char > &I)
vpTemplateTrackerWarp * Warp
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:5071
void computeProba(int &nbpoint)
virtual void dWarp(const vpColVector &X1, const vpColVector &X2, const vpColVector &ParamM, vpMatrix &dW)=0
vpHessienType hessianComputation