Visual Servoing Platform  version 3.2.1 under development (2019-08-21) under development (2019-08-21)
vpTemplateTrackerSSDForwardAdditional.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  * Template tracker.
33  *
34  * Authors:
35  * Amaury Dame
36  * Aurelien Yol
37  * Fabien Spindler
38  *
39  *****************************************************************************/
40 
41 #include <limits> // numeric_limits
42 
43 #include <visp3/core/vpImageTools.h>
44 #include <visp3/tt/vpTemplateTrackerSSDForwardAdditional.h>
45 
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  do {
69  unsigned int Nbpoint = 0;
70  double erreur = 0;
71  G = 0;
72  H = 0;
73  Warp->computeCoeff(p);
74  for (unsigned int point = 0; point < templateSize; point++) {
75  i = ptTemplate[point].y;
76  j = ptTemplate[point].x;
77  X1[0] = j;
78  X1[1] = i;
79 
80  Warp->computeDenom(X1, p);
81  Warp->warpX(X1, X2, p);
82 
83  j2 = X2[0];
84  i2 = X2[1];
85  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
86  Tij = ptTemplate[point].val;
87 
88  if (!blur)
89  IW = I.getValue(i2, j2);
90  else
91  IW = BI.getValue(i2, j2);
92 
93  dIWx = dIx.getValue(i2, j2);
94  dIWy = dIy.getValue(i2, j2);
95  Nbpoint++;
96  // Calcul du Hessien
97  Warp->dWarp(X1, X2, p, dW);
98  double *tempt = new double[nbParam];
99  for (unsigned int it = 0; it < nbParam; it++)
100  tempt[it] = dW[0][it] * dIWx + dW[1][it] * dIWy;
101 
102  for (unsigned int it = 0; it < nbParam; it++)
103  for (unsigned int jt = 0; jt < nbParam; jt++)
104  H[it][jt] += tempt[it] * tempt[jt];
105 
106  double er = (Tij - IW);
107  for (unsigned int it = 0; it < nbParam; it++)
108  G[it] += er * tempt[it];
109 
110  erreur += (er * er);
111  delete[] tempt;
112  }
113  }
114  if (Nbpoint == 0) {
115  // std::cout<<"plus de point dans template suivi"<<std::endl;
116  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
117  }
118 
119  vpMatrix::computeHLM(H, lambda, HLM);
120  try {
121  dp = 1. * HLM.inverseByLU() * G;
122  } catch (const vpException &e) {
123  throw(e);
124  }
125 
126  switch (minimizationMethod) {
128  vpColVector p_test_LMA(nbParam);
129  p_test_LMA = p + 1. * dp;
130  erreur = -getCost(I, p);
131  double erreur_LMA = -getCost(I, p_test_LMA);
132  if (erreur_LMA < erreur) {
133  p = p_test_LMA;
134  lambda = (lambda / 10. < 1e-6) ? lambda / 10. : 1e-6;
135  } else {
136  lambda = (lambda * 10. < 1e6) ? 1e6 : lambda * 10.;
137  }
138  } break;
140  dp = gain * 0.000001 * G;
141  if (useBrent) {
142  alpha = 2.;
143  computeOptimalBrentGain(I, p, erreur, dp, alpha);
144  dp = alpha * dp;
145  }
146  p += 1. * dp;
147  break;
148  }
149 
151  if (iterationGlobale != 0) {
152  vpColVector s_quasi = p - p_prec;
153  vpColVector y_quasi = G - G_prec;
154  double s_scal_y = s_quasi.t() * y_quasi;
155  // if(s_scal_y!=0)//BFGS
156  // 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;
157  // if(s_scal_y!=0.0)//DFP
158  if (std::fabs(s_scal_y) > std::numeric_limits<double>::epsilon()) // DFP
159  KQuasiNewton = KQuasiNewton + 0.001 * (s_quasi * s_quasi.t() / s_scal_y -
160  KQuasiNewton * y_quasi * y_quasi.t() * KQuasiNewton /
161  (y_quasi.t() * KQuasiNewton * y_quasi));
162  }
163  dp = -KQuasiNewton * G;
164  p_prec = p;
165  G_prec = G;
166  p -= 1.01 * dp;
167  } break;
168 
170  default: {
171  if (useBrent) {
172  alpha = 2.;
173  computeOptimalBrentGain(I, p, erreur, dp, alpha);
174  dp = alpha * dp;
175  }
176 
177  p += 1. * dp;
178  break;
179  }
180  }
181 
182  iteration++;
184  } while (/*( erreur_prec-erreur<50) && */ (iteration < iterationMax));
185 
186  // std::cout<<"erreur "<<erreur<<std::endl;
187  nbIteration = iteration;
188 }
void trackNoPyr(const vpImage< unsigned char > &I)
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1422
vpTemplateTrackerPoint * ptTemplate
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
vpMatrix inverseByLU() const
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
vpRowVector t() 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
unsigned int templateSize
unsigned int iterationMax
double getCost(const vpImage< unsigned char > &I, const vpColVector &tp)
Error that can be emited by the vpTracker class and its derivates.
vpImage< double > dIx
vpImage< double > dIy
unsigned int iterationGlobale
unsigned int nbIteration
unsigned int getHeight() const
Definition: vpImage.h:186
Implementation of column vector and the associated operations.
Definition: vpColVector.h:130
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, const bool convolve=false)
vpTemplateTrackerWarp * Warp
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:5233
unsigned int getWidth() const
Definition: vpImage.h:244
virtual void dWarp(const vpColVector &X1, const vpColVector &X2, const vpColVector &ParamM, vpMatrix &dW)=0