Visual Servoing Platform  version 3.6.1 under development (2024-10-10)
vpTemplateTrackerMIForwardAdditional.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2023 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 https://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  *
38 *****************************************************************************/
39 
40 #include <visp3/tt_mi/vpTemplateTrackerMIForwardAdditional.h>
41 
42 #ifdef VISP_HAVE_OPENMP
43 #include <omp.h>
44 #endif
45 
46 BEGIN_VISP_NAMESPACE
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  }
127  catch (const vpException &e) {
128  throw(e);
129  }
130  }
131 }
132 
134 {
135  dW = 0;
136 
137  int Nbpoint = 0;
138  if (blur)
142 
143  double MI = 0, MIprec = -1000;
144 
145  MI_preEstimation = -getCost(I, p);
146 
147  double alpha = 2.;
148 
149  unsigned int iteration = 0;
150 
151  initPosEvalRMS(p);
152  double evolRMS_init = 0;
153  double evolRMS_prec = 0;
154  double evolRMS_delta;
155 
156  do {
157  if (iteration % 5 == 0)
159  Nbpoint = 0;
160  MIprec = MI;
161  MI = 0;
162  // erreur=0;
163 
165 
166  Warp->computeCoeff(p);
167 #ifdef VISP_HAVE_OPENMP
168  int nthreads = omp_get_num_procs();
169  // std::cout << "file: " __FILE__ << " line: " << __LINE__ << " function:
170  // " << __FUNCTION__ << " nthread: " << nthreads << std::endl;
171  omp_set_num_threads(nthreads);
172 #pragma omp parallel for default(shared)
173 #endif
174  for (int point = 0; point < (int)templateSize; point++) {
175  int i = ptTemplate[point].y;
176  int j = ptTemplate[point].x;
177  X1[0] = j;
178  X1[1] = i;
179 
180  Warp->computeDenom(X1, p);
181  Warp->warpX(X1, X2, p);
182 
183  double j2 = X2[0];
184  double i2 = X2[1];
185 
186  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
187  Nbpoint++;
188  double Tij = ptTemplate[point].val;
189  double IW;
190  if (!blur)
191  IW = I.getValue(i2, j2);
192  else
193  IW = BI.getValue(i2, j2);
194 
195  double dx = dIx.getValue(i2, j2) * (Nc - 1) / 255.;
196  double dy = dIy.getValue(i2, j2) * (Nc - 1) / 255.;
197 
198  int ct = (int)((IW * (Nc - 1)) / 255.);
199  int cr = (int)((Tij * (Nc - 1)) / 255.);
200  double et = (IW * (Nc - 1)) / 255. - ct;
201  double er = ((double)Tij * (Nc - 1)) / 255. - cr;
202 
203  Warp->dWarp(X1, X2, p, dW);
204 
205  double *tptemp = new double[nbParam];
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  }
222  else {
223  computeProba(Nbpoint);
224  computeMI(MI);
225  computeHessien(H);
226  computeGradient();
227 
229  try {
230  switch (hessianComputation) {
232  dp = gain * HLMdesireInverse * G;
233  break;
235  if (HLM.cond() > HLMdesire.cond())
236  dp = gain * HLMdesireInverse * G;
237  else
238  dp = gain * 0.2 * HLM.inverseByLU() * G;
239  break;
240  default:
241  dp = gain * 0.2 * HLM.inverseByLU() * G;
242  break;
243  }
244  }
245  catch (const vpException &e) {
246  throw(e);
247  }
248  }
249 
250  switch (minimizationMethod) {
252  vpColVector p_test_LMA(nbParam);
254  p_test_LMA = p - 100000.1 * dp;
255  else
256  p_test_LMA = p + dp;
257  MI = -getCost(I, p);
258  double MI_LMA = -getCost(I, p_test_LMA);
259  if (MI_LMA > MI) {
260  p = p_test_LMA;
261  lambda = (lambda / 10. < 1e-6) ? lambda / 10. : 1e-6;
262  }
263  else {
264  lambda = (lambda * 10. < 1e6) ? 1e6 : lambda * 10.;
265  }
266  } break;
268  dp = -gain * 6.0 * G;
269  if (useBrent) {
270  alpha = 2.;
271  computeOptimalBrentGain(I, p, -MI, dp, alpha);
272  dp = alpha * dp;
273  }
274  p += dp;
275  break;
276  }
277 
279  if (iterationGlobale != 0) {
280  vpColVector s_quasi = p - p_prec;
281  vpColVector y_quasi = G - G_prec;
282  double s_scal_y = s_quasi.t() * y_quasi;
283  if (std::fabs(s_scal_y) > std::numeric_limits<double>::epsilon()) {
284  KQuasiNewton = KQuasiNewton + 0.001 * (s_quasi * s_quasi.t() / s_scal_y -
285  KQuasiNewton * y_quasi * y_quasi.t() * KQuasiNewton /
286  (y_quasi.t() * KQuasiNewton * y_quasi));
287  }
288  }
289  dp = -KQuasiNewton * G;
290  p_prec = p;
291  G_prec = G;
292  p -= 1.01 * dp;
293  } break;
294 
295  default: {
297  dp = -0.1 * dp;
298  if (useBrent) {
299  alpha = 2.;
300  computeOptimalBrentGain(I, p, -MI, dp, alpha);
301  dp = alpha * dp;
302  }
303 
304  p += dp;
305  break;
306  }
307  }
308 
309  computeEvalRMS(p);
310 
311  if (iteration == 0) {
312  evolRMS_init = evolRMS;
313  }
314  iteration++;
316 
317  evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
318  evolRMS_prec = evolRMS;
319 
320  } while ((std::fabs(MI - MIprec) > std::fabs(MI) * std::numeric_limits<double>::epsilon()) &&
321  (iteration < iterationMax) && (evolRMS_delta > std::fabs(evolRMS_init) * evolRMS_eps));
322 
323  if (Nbpoint == 0) {
324  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
325  }
326 
327  nbIteration = iteration;
328  MI_postEstimation = -getCost(I, p);
330  MI_postEstimation = -1;
331  }
332 }
333 END_VISP_NAMESPACE
Implementation of column vector and the associated operations.
Definition: vpColVector.h:191
vpRowVector t() const
error that can be emitted by ViSP classes.
Definition: vpException.h:60
static void getGradXGauss2D(const vpImage< ImageType > &I, vpImage< FilterType > &dIx, const FilterType *gaussianKernel, const FilterType *gaussianDerivativeKernel, unsigned int size, const vpImage< bool > *p_mask=nullptr)
static void filter(const vpImage< ImageType > &I, vpImage< FilterType > &If, const vpArray2D< FilterType > &M, bool convolve=false, const vpImage< bool > *p_mask=nullptr)
static void getGradYGauss2D(const vpImage< ImageType > &I, vpImage< FilterType > &dIy, const FilterType *gaussianKernel, const FilterType *gaussianDerivativeKernel, unsigned int size, const vpImage< bool > *p_mask=nullptr)
unsigned int getWidth() const
Definition: vpImage.h:242
Type getValue(unsigned int i, unsigned int j) const
unsigned int getHeight() const
Definition: vpImage.h:181
double cond(double svThreshold=1e-6) const
Definition: vpMatrix.cpp:1810
vpMatrix inverseByLU() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:1873
void initHessienDesired(const vpImage< unsigned char > &I)
void trackNoPyr(const vpImage< unsigned char > &I)
vpHessienApproximationType ApproxHessian
vpHessienType hessianComputation
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp) VP_OVERRIDE
void computeMI(double &MI)
void computeProba(int &nbpoint)
void computeHessien(vpMatrix &H)
virtual void dWarp(const vpColVector &X1, const vpColVector &X2, const vpColVector &p, vpMatrix &dM)=0
virtual void warpX(const int &v1, const int &u1, double &v2, double &u2, const vpColVector &p)=0
vpImage< double > dIx
vpImage< double > dIy
void computeEvalRMS(const vpColVector &p)
void computeOptimalBrentGain(const vpImage< unsigned char > &I, vpColVector &tp, double tMI, vpColVector &direction, double &alpha)
unsigned int iterationMax
void initPosEvalRMS(const vpColVector &p)
unsigned int nbIteration
vpTemplateTrackerPoint * ptTemplate
vpTemplateTrackerWarp * Warp
unsigned int iterationGlobale
vpImage< double > BI
unsigned int templateSize
Error that can be emitted by the vpTracker class and its derivatives.
@ notEnoughPointError
Not enough point to track.