Visual Servoing Platform  version 3.6.1 under development (2024-12-12)
vpTemplateTrackerSSDForwardAdditional.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  * Template tracker.
33  *
34  * Authors:
35  * Amaury Dame
36  * Aurelien Yol
37  *
38 *****************************************************************************/
39 
40 #include <limits> // numeric_limits
41 
42 #include <visp3/core/vpImageTools.h>
43 #include <visp3/tt/vpTemplateTrackerSSDForwardAdditional.h>
44 
45 BEGIN_VISP_NAMESPACE
47  : vpTemplateTrackerSSD(warp), minimizationMethod(USE_NEWTON), p_prec(), G_prec(), KQuasiNewton()
48 {
49  useCompositionnal = false;
50 }
51 
53 {
54  if (blur)
58 
59  dW = 0;
60 
61  double lambda = lambdaDep;
62  double IW, dIWx, dIWy;
63  double Tij;
64  unsigned int iteration = 0;
65  int i, j;
66  double i2, j2;
67  double alpha = 2.;
68 
70 
71  double evolRMS_init = 0;
72  double evolRMS_prec = 0;
73  double evolRMS_delta;
74  double *tempt = new double[nbParam];
75 
76  do {
77  unsigned int Nbpoint = 0;
78  double erreur = 0;
79  G = 0;
80  H = 0;
81  Warp->computeCoeff(p);
82  for (unsigned int point = 0; point < templateSize; point++) {
83  i = ptTemplate[point].y;
84  j = ptTemplate[point].x;
85  X1[0] = j;
86  X1[1] = i;
87 
88  Warp->computeDenom(X1, p);
89  Warp->warpX(X1, X2, p);
90 
91  j2 = X2[0];
92  i2 = X2[1];
93  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
94  Tij = ptTemplate[point].val;
95 
96  if (!blur)
97  IW = I.getValue(i2, j2);
98  else
99  IW = BI.getValue(i2, j2);
100 
101  dIWx = dIx.getValue(i2, j2);
102  dIWy = dIy.getValue(i2, j2);
103  Nbpoint++;
104  // Calcul du Hessien
105  Warp->dWarp(X1, X2, p, dW);
106  for (unsigned int it = 0; it < nbParam; it++)
107  tempt[it] = dW[0][it] * dIWx + dW[1][it] * dIWy;
108 
109  for (unsigned int it = 0; it < nbParam; it++)
110  for (unsigned int jt = 0; jt < nbParam; jt++)
111  H[it][jt] += tempt[it] * tempt[jt];
112 
113  double er = (Tij - IW);
114  for (unsigned int it = 0; it < nbParam; it++)
115  G[it] += er * tempt[it];
116 
117  erreur += (er * er);
118  }
119  }
120  if (Nbpoint == 0) {
121  delete[] tempt;
122  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
123  }
124 
125  vpMatrix::computeHLM(H, lambda, HLM);
126  try {
127  dp = HLM.inverseByLU() * G;
128  }
129  catch (const vpException &e) {
130  delete[] tempt;
131  throw(e);
132  }
133 
134  switch (minimizationMethod) {
136  vpColVector p_test_LMA(nbParam);
137  p_test_LMA = p + dp;
138  erreur = -getCost(I, p);
139  double erreur_LMA = -getCost(I, p_test_LMA);
140  if (erreur_LMA < erreur) {
141  p = p_test_LMA;
142  lambda = (lambda / 10. < 1e-6) ? lambda / 10. : 1e-6;
143  }
144  else {
145  lambda = (lambda * 10. < 1e6) ? 1e6 : lambda * 10.;
146  }
147  } break;
149  dp = gain * 0.000001 * G;
150  if (useBrent) {
151  alpha = 2.;
152  computeOptimalBrentGain(I, p, erreur, dp, alpha);
153  dp = alpha * dp;
154  }
155  p += dp;
156  break;
157  }
158 
160  if (iterationGlobale != 0) {
161  vpColVector s_quasi = p - p_prec;
162  vpColVector y_quasi = G - G_prec;
163  double s_scal_y = s_quasi.t() * y_quasi;
164  if (std::fabs(s_scal_y) > std::numeric_limits<double>::epsilon()) // DFP
165  KQuasiNewton = KQuasiNewton + 0.001 * (s_quasi * s_quasi.t() / s_scal_y -
166  KQuasiNewton * y_quasi * y_quasi.t() * KQuasiNewton /
167  (y_quasi.t() * KQuasiNewton * y_quasi));
168  }
169  dp = -KQuasiNewton * G;
170  p_prec = p;
171  G_prec = G;
172  p -= 1.01 * dp;
173  } break;
174 
176  default: {
177  if (useBrent) {
178  alpha = 2.;
179  computeOptimalBrentGain(I, p, erreur, dp, alpha);
180  dp = alpha * dp;
181  }
182 
183  p += dp;
184  break;
185  }
186  }
187 
188  computeEvalRMS(p);
189 
190  if (iteration == 0) {
191  evolRMS_init = evolRMS;
192  }
193 
194  iteration++;
196 
197  evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
198  evolRMS_prec = evolRMS;
199 
200  } while ((iteration < iterationMax) && (evolRMS_delta > std::fabs(evolRMS_init) * evolRMS_eps));
201  delete[] tempt;
202 
203  nbIteration = iteration;
204 }
205 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
vpMatrix inverseByLU() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:1869
void trackNoPyr(const vpImage< unsigned char > &I)
VP_EXPLICIT vpTemplateTrackerSSDForwardAdditional(vpTemplateTrackerWarp *warp)
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp)
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.