Visual Servoing Platform  version 3.3.0
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), p_prec(), G_prec(), KQuasiNewton()
49 {
50  useCompositionnal = false;
51 }
52 
54 {
55  dW = 0;
56 
57  int Nbpoint = 0;
58 
59  if (blur)
63 
64  double Tij;
65  double IW, dx, dy;
66  int cr, ct;
67  double er, et;
68 
69  Nbpoint = 0;
70 
72  Warp->computeCoeff(p);
73  for (unsigned int point = 0; point < templateSize; point++) {
74  int i = ptTemplate[point].y;
75  int j = ptTemplate[point].x;
76  X1[0] = j;
77  X1[1] = i;
78  X2[0] = j;
79  X2[1] = i;
80 
81  Warp->computeDenom(X1, p);
82  Warp->warpX(X1, X2, p);
83 
84  double j2 = X2[0];
85  double i2 = X2[1];
86 
87  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
88  Nbpoint++;
89  Tij = ptTemplate[point].val;
90  if (!blur)
91  IW = I.getValue(i2, j2);
92  else
93  IW = BI.getValue(i2, j2);
94 
95  dx = dIx.getValue(i2, j2) * (Nc - 1) / 255.;
96  dy = dIy.getValue(i2, j2) * (Nc - 1) / 255.;
97 
98  ct = static_cast<int>((IW * (Nc - 1)) / 255.);
99  cr = static_cast<int>((Tij * (Nc - 1)) / 255.);
100  et = (IW * (Nc - 1)) / 255. - ct;
101  er = (static_cast<double>(Tij) * (Nc - 1)) / 255. - cr;
102  Warp->dWarp(X1, X2, p, dW);
103 
104  double *tptemp = new double[nbParam];
105  for (unsigned int it = 0; it < nbParam; it++)
106  tptemp[it] = dW[0][it] * dx + dW[1][it] * dy;
107 
109  vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
111  vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
112 
113  delete[] tptemp;
114  }
115  }
116 
117  if (Nbpoint > 0) {
118  double MI;
119  computeProba(Nbpoint);
120  computeMI(MI);
122 
124  try {
126  } catch (const vpException &e) {
127  throw(e);
128  }
129  }
130 }
131 
133 {
134  dW = 0;
135 
136  int Nbpoint = 0;
137  if (blur)
141 
142  double MI = 0, MIprec = -1000;
143 
144  MI_preEstimation = -getCost(I, p);
145 
146  double alpha = 2.;
147 
148  unsigned int iteration = 0;
149 
150  initPosEvalRMS(p);
151  double evolRMS_init = 0;
152  double evolRMS_prec = 0;
153  double evolRMS_delta;
154 
155  do {
156  if (iteration % 5 == 0)
158  Nbpoint = 0;
159  MIprec = MI;
160  MI = 0;
161  // erreur=0;
162 
164 
165  Warp->computeCoeff(p);
166 #ifdef VISP_HAVE_OPENMP
167  int nthreads = omp_get_num_procs();
168  // std::cout << "file: " __FILE__ << " line: " << __LINE__ << " function:
169  // " << __FUNCTION__ << " nthread: " << nthreads << std::endl;
170  omp_set_num_threads(nthreads);
171 #pragma omp parallel for default(shared)
172 #endif
173  for (int point = 0; point < (int)templateSize; point++) {
174  int i = ptTemplate[point].y;
175  int j = ptTemplate[point].x;
176  X1[0] = j;
177  X1[1] = i;
178 
179  Warp->computeDenom(X1, p);
180  Warp->warpX(X1, X2, p);
181 
182  double j2 = X2[0];
183  double i2 = X2[1];
184 
185  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
186  Nbpoint++;
187  double Tij = ptTemplate[point].val;
188  double IW;
189  if (!blur)
190  IW = I.getValue(i2, j2);
191  else
192  IW = BI.getValue(i2, j2);
193 
194  double dx = dIx.getValue(i2, j2) * (Nc - 1) / 255.;
195  double dy = dIy.getValue(i2, j2) * (Nc - 1) / 255.;
196 
197  int ct = (int)((IW * (Nc - 1)) / 255.);
198  int cr = (int)((Tij * (Nc - 1)) / 255.);
199  double et = (IW * (Nc - 1)) / 255. - ct;
200  double er = ((double)Tij * (Nc - 1)) / 255. - cr;
201 
202  Warp->dWarp(X1, X2, p, dW);
203 
204  double *tptemp = new double[nbParam];
205  ;
206  for (unsigned int it = 0; it < nbParam; it++)
207  tptemp[it] = (dW[0][it] * dx + dW[1][it] * dy);
209  vpTemplateTrackerMIBSpline::PutTotPVBsplineNoSecond(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
211  vpTemplateTrackerMIBSpline::PutTotPVBspline(PrtTout, cr, er, ct, et, Nc, tptemp, nbParam, bspline);
212 
213  delete[] tptemp;
214  }
215  }
216 
217  if (Nbpoint == 0) {
218  diverge = true;
219  MI = 0;
220  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
221  } else {
222  computeProba(Nbpoint);
223  computeMI(MI);
224  computeHessien(H);
225  computeGradient();
226 
228  try {
229  switch (hessianComputation) {
231  dp = gain * HLMdesireInverse * G;
232  break;
234  if (HLM.cond() > HLMdesire.cond())
235  dp = gain * HLMdesireInverse * G;
236  else
237  dp = gain * 0.2 * HLM.inverseByLU() * G;
238  break;
239  default:
240  dp = gain * 0.2 * HLM.inverseByLU() * G;
241  break;
242  }
243  } catch (const vpException &e) {
244  throw(e);
245  }
246  }
247 
248  switch (minimizationMethod) {
250  vpColVector p_test_LMA(nbParam);
252  p_test_LMA = p - 100000.1 * dp;
253  else
254  p_test_LMA = p + dp;
255  MI = -getCost(I, p);
256  double MI_LMA = -getCost(I, p_test_LMA);
257  if (MI_LMA > MI) {
258  p = p_test_LMA;
259  lambda = (lambda / 10. < 1e-6) ? lambda / 10. : 1e-6;
260  } else {
261  lambda = (lambda * 10. < 1e6) ? 1e6 : lambda * 10.;
262  }
263  } break;
265  dp = -gain * 6.0 * G;
266  if (useBrent) {
267  alpha = 2.;
268  computeOptimalBrentGain(I, p, -MI, dp, alpha);
269  dp = alpha * dp;
270  }
271  p += dp;
272  break;
273  }
274 
276  if (iterationGlobale != 0) {
277  vpColVector s_quasi = p - p_prec;
278  vpColVector y_quasi = G - G_prec;
279  double s_scal_y = s_quasi.t() * y_quasi;
280  if (std::fabs(s_scal_y) > std::numeric_limits<double>::epsilon()) {
281  KQuasiNewton = KQuasiNewton + 0.001 * (s_quasi * s_quasi.t() / s_scal_y -
282  KQuasiNewton * y_quasi * y_quasi.t() * KQuasiNewton /
283  (y_quasi.t() * KQuasiNewton * y_quasi));
284  }
285  }
286  dp = -KQuasiNewton * G;
287  p_prec = p;
288  G_prec = G;
289  p -= 1.01 * dp;
290  } break;
291 
292  default: {
294  dp = -0.1 * dp;
295  if (useBrent) {
296  alpha = 2.;
297  computeOptimalBrentGain(I, p, -MI, dp, alpha);
298  dp = alpha * dp;
299  }
300 
301  p += dp;
302  break;
303  }
304  }
305 
306  computeEvalRMS(p);
307 
308  if (iteration == 0) {
309  evolRMS_init = evolRMS;
310  }
311  iteration++;
313 
314  evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
315  evolRMS_prec = evolRMS;
316 
317  } while ((std::fabs(MI - MIprec) > std::fabs(MI) * std::numeric_limits<double>::epsilon()) &&
318  (iteration < iterationMax) && (evolRMS_delta > std::fabs(evolRMS_init)*evolRMS_eps));
319 
320  if (Nbpoint == 0) {
321  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
322  }
323 
324  nbIteration = iteration;
325  MI_postEstimation = -getCost(I, p);
327  MI_postEstimation = -1;
328  }
329 }
vpTrackingException
Error that can be emited by the vpTracker class and its derivates.
Definition: vpTrackingException.h:68
vpTemplateTracker::G
vpColVector G
Definition: vpTemplateTracker.h:104
vpTemplateTracker::taillef
unsigned int taillef
Definition: vpTemplateTracker.h:112
vpImage::getValue
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1457
vpTemplateTrackerWarp
Definition: vpTemplateTrackerWarp.h:57
vpTemplateTrackerMIForwardAdditional::initHessienDesired
void initHessienDesired(const vpImage< unsigned char > &I)
Definition: vpTemplateTrackerMIForwardAdditional.cpp:52
vpTemplateTracker::initPosEvalRMS
void initPosEvalRMS(const vpColVector &p)
Definition: vpTemplateTracker.cpp:957
vpTemplateTracker::ptTemplate
vpTemplateTrackerPoint * ptTemplate
Definition: vpTemplateTracker.h:74
vpColVector::t
vpRowVector t() const
Definition: vpColVector.cpp:781
vpTemplateTracker::HLM
vpMatrix HLM
Definition: vpTemplateTracker.h:99
vpTemplateTrackerMI::computeGradient
void computeGradient()
Definition: vpTemplateTrackerMI.cpp:556
vpTemplateTracker::useBrent
bool useBrent
Definition: vpTemplateTracker.h:110
vpTemplateTracker::dp
vpColVector dp
Definition: vpTemplateTracker.h:132
vpTemplateTrackerMIForwardAdditional::USE_QUASINEWTON
Definition: vpTemplateTrackerMIForwardAdditional.h:61
vpTemplateTracker::diverge
bool diverge
Definition: vpTemplateTracker.h:124
vpImageFilter::filter
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, bool convolve=false)
Definition: vpImageFilter.cpp:126
vpTemplateTrackerMIForwardAdditional::USE_LMA
Definition: vpTemplateTrackerMIForwardAdditional.h:61
vpTemplateTracker::X2
vpColVector X2
Definition: vpTemplateTracker.h:136
vpTemplateTrackerMI::computeProba
void computeProba(int &nbpoint)
Definition: vpTemplateTrackerMI.cpp:342
vpTemplateTracker::fgdG
double * fgdG
Definition: vpTemplateTracker.h:114
vpTemplateTracker::useCompositionnal
bool useCompositionnal
Definition: vpTemplateTracker.h:126
vpTemplateTracker::computeOptimalBrentGain
void computeOptimalBrentGain(const vpImage< unsigned char > &I, vpColVector &tp, double tMI, vpColVector &direction, double &alpha)
Definition: vpTemplateTracker.cpp:391
vpTemplateTracker::fgG
double * fgG
Definition: vpTemplateTracker.h:113
vpTemplateTrackerMI
Definition: vpTemplateTrackerMI.h:52
vpTemplateTrackerMI::HESSIAN_0
Definition: vpTemplateTrackerMI.h:58
vpTemplateTracker::evolRMS
double evolRMS
Definition: vpTemplateTracker.h:69
vpTemplateTrackerMIForwardAdditional::trackNoPyr
void trackNoPyr(const vpImage< unsigned char > &I)
Definition: vpTemplateTrackerMIForwardAdditional.cpp:131
vpTemplateTracker::blur
bool blur
Definition: vpTemplateTracker.h:109
vpTemplateTracker::HLMdesire
vpMatrix HLMdesire
Definition: vpTemplateTracker.h:100
vpTemplateTrackerMI::MI_postEstimation
double MI_postEstimation
Definition: vpTemplateTrackerMI.h:99
vpTrackingException::notEnoughPointError
Definition: vpTrackingException.h:79
vpTemplateTrackerWarp::dWarp
virtual void dWarp(const vpColVector &X1, const vpColVector &X2, const vpColVector &ParamM, vpMatrix &dW)=0
vpTemplateTrackerMI::computeMI
void computeMI(double &MI)
Definition: vpTemplateTrackerMI.cpp:387
vpTemplateTrackerMI::USE_HESSIEN_BEST_COND
Definition: vpTemplateTrackerMI.h:66
vpImageFilter::getGradYGauss2D
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
Definition: vpImageFilter.cpp:766
vpTemplateTrackerWarp::warpX
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
vpTemplateTracker::templateSize
unsigned int templateSize
Definition: vpTemplateTracker.h:77
vpTemplateTrackerMI::HESSIAN_NEW
Definition: vpTemplateTrackerMI.h:62
vpMatrix::cond
double cond(double svThreshold=1e-6) const
Definition: vpMatrix.cpp:5391
vpTemplateTracker::computeEvalRMS
void computeEvalRMS(const vpColVector &p)
Definition: vpTemplateTracker.cpp:924
vpTemplateTracker::X1
vpColVector X1
Definition: vpTemplateTracker.h:135
vpColVector
Implementation of column vector and the associated operations.
Definition: vpColVector.h:129
vpTemplateTrackerMI::MI_preEstimation
double MI_preEstimation
Definition: vpTemplateTrackerMI.h:98
vpTemplateTracker::BI
vpImage< double > BI
Definition: vpTemplateTracker.h:140
vpTemplateTracker::nbIteration
unsigned int nbIteration
Definition: vpTemplateTracker.h:125
vpTemplateTracker::dW
vpMatrix dW
Definition: vpTemplateTracker.h:138
vpTemplateTracker::Hdesire
vpMatrix Hdesire
Definition: vpTemplateTracker.h:97
vpTemplateTrackerPoint::val
double val
Definition: vpTemplateTrackerHeader.h:74
vpTemplateTracker::iterationMax
unsigned int iterationMax
Definition: vpTemplateTracker.h:120
vpTemplateTracker::HLMdesireInverse
vpMatrix HLMdesireInverse
Definition: vpTemplateTracker.h:102
vpTemplateTrackerMI::PrtTout
double * PrtTout
Definition: vpTemplateTrackerMI.h:82
vpTemplateTracker::nbParam
unsigned int nbParam
Definition: vpTemplateTracker.h:118
vpImage::getHeight
unsigned int getHeight() const
Definition: vpImage.h:221
vpTemplateTracker::dIy
vpImage< double > dIy
Definition: vpTemplateTracker.h:142
vpTemplateTrackerMIForwardAdditional::vpTemplateTrackerMIForwardAdditional
vpTemplateTrackerMIForwardAdditional()
Default constructor.
Definition: vpTemplateTrackerMIForwardAdditional.h:76
vpMatrix::inverseByLU
vpMatrix inverseByLU() const
Definition: vpMatrix_lu.cpp:133
vpTemplateTracker::H
vpMatrix H
Definition: vpTemplateTracker.h:96
vpTemplateTrackerMI::ApproxHessian
vpHessienApproximationType ApproxHessian
Definition: vpTemplateTrackerMI.h:73
vpTemplateTracker::iterationGlobale
unsigned int iterationGlobale
Definition: vpTemplateTracker.h:122
vpTemplateTrackerMI::USE_HESSIEN_DESIRE
Definition: vpTemplateTrackerMI.h:66
vpTemplateTrackerPoint::y
int y
Definition: vpTemplateTrackerHeader.h:72
vpTemplateTracker::dIx
vpImage< double > dIx
Definition: vpTemplateTracker.h:141
vpTemplateTrackerMI::bspline
int bspline
Definition: vpTemplateTrackerMI.h:89
vpTemplateTrackerMI::zeroProbabilities
void zeroProbabilities()
Definition: vpTemplateTrackerMI.cpp:584
vpTemplateTrackerMIForwardAdditional::USE_GRADIENT
Definition: vpTemplateTrackerMIForwardAdditional.h:61
vpTemplateTrackerMI::getCost
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp)
Definition: vpTemplateTrackerMI.cpp:147
vpTemplateTracker::p
vpColVector p
Definition: vpTemplateTracker.h:131
vpTemplateTrackerMI::HESSIAN_NONSECOND
Definition: vpTemplateTrackerMI.h:57
vpMatrix::computeHLM
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:5452
vpTemplateTrackerMI::hessianComputation
vpHessienType hessianComputation
Definition: vpTemplateTrackerMI.h:72
vpTemplateTrackerMI::Nc
int Nc
Definition: vpTemplateTrackerMI.h:91
vpImage< unsigned char >
vpTemplateTracker::gain
double gain
Definition: vpTemplateTracker.h:106
vpTemplateTrackerPoint::x
int x
Definition: vpTemplateTrackerHeader.h:72
vpTemplateTrackerMI::lambda
double lambda
Definition: vpTemplateTrackerMI.h:74
vpException
error that can be emited by ViSP classes.
Definition: vpException.h:70
vpTemplateTrackerMI::computeHessien
void computeHessien(vpMatrix &H)
Definition: vpTemplateTrackerMI.cpp:420
vpTemplateTracker::Warp
vpTemplateTrackerWarp * Warp
Definition: vpTemplateTracker.h:129
vpImage::getWidth
unsigned int getWidth() const
Definition: vpImage.h:279
vpImageFilter::getGradXGauss2D
static void getGradXGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIx, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
Definition: vpImageFilter.cpp:749
vpTemplateTracker::evolRMS_eps
double evolRMS_eps
Definition: vpTemplateTracker.h:72