Visual Servoing Platform  version 3.6.1 under development (2024-10-15)
vpTemplateTrackerSSDESM.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 #include <visp3/core/vpImageFilter.h>
40 #include <visp3/tt/vpTemplateTrackerSSDESM.h>
41 
42 BEGIN_VISP_NAMESPACE
44  : vpTemplateTrackerSSD(warp), compoInitialised(false), HDir(), HInv(), HLMDir(), HLMInv(), GDir(), GInv()
45 {
46  useCompositionnal = false;
47  useInverse = false;
48 
49  if (!Warp->isESMcompatible()) {
50  throw(vpException(vpException::badValue, "The selected warp function is not appropriate for the ESM algorithm..."));
51  }
52 
59 }
60 
62 
64 {
66  int i, j;
67  // direct
68  for (unsigned int point = 0; point < templateSize; point++) {
69  i = ptTemplate[point].y;
70  j = ptTemplate[point].x;
71  ptTemplateCompo[point].dW = new double[2 * nbParam];
72  Warp->getdWdp0(i, j, ptTemplateCompo[point].dW);
73  }
74 
75  // inverse
76  HInv = 0;
77  for (unsigned int point = 0; point < templateSize; point++) {
78  i = ptTemplate[point].y;
79  j = ptTemplate[point].x;
80 
81  ptTemplate[point].dW = new double[nbParam];
82  Warp->getdW0(i, j, ptTemplate[point].dy, ptTemplate[point].dx, ptTemplate[point].dW);
83 
84  for (unsigned int it = 0; it < nbParam; it++)
85  for (unsigned int jt = 0; jt < nbParam; jt++)
86  HInv[it][jt] += ptTemplate[point].dW[it] * ptTemplate[point].dW[jt];
87  }
89 
90  compoInitialised = true;
91 }
92 
94 {
95  if (blur)
99 
100  double IW, dIWx, dIWy;
101  double Tij;
102  unsigned int iteration = 0;
103  int i, j;
104  double i2, j2;
105  double alpha = 2.;
106 
107  initPosEvalRMS(p);
108 
109  double evolRMS_init = 0;
110  double evolRMS_prec = 0;
111  double evolRMS_delta;
112  double *tempt = new double[nbParam];
113 
114  do {
115  unsigned int Nbpoint = 0;
116  double erreur = 0;
117  dp = 0;
118  HDir = 0;
119  GDir = 0;
120  GInv = 0;
121  Warp->computeCoeff(p);
122  for (unsigned int point = 0; point < templateSize; point++) {
123  i = ptTemplate[point].y;
124  j = ptTemplate[point].x;
125  X1[0] = j;
126  X1[1] = i;
127 
128  Warp->computeDenom(X1, p);
129  Warp->warpX(X1, X2, p);
130 
131  j2 = X2[0];
132  i2 = X2[1];
133  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
134  // INVERSE
135  Tij = ptTemplate[point].val;
136  if (!blur)
137  IW = I.getValue(i2, j2);
138  else
139  IW = BI.getValue(i2, j2);
140  Nbpoint++;
141  double er = (Tij - IW);
142  for (unsigned int it = 0; it < nbParam; it++)
143  GInv[it] += er * ptTemplate[point].dW[it];
144 
145  erreur += er * er;
146 
147  dIWx = dIx.getValue(i2, j2) + ptTemplate[point].dx;
148  dIWy = dIy.getValue(i2, j2) + ptTemplate[point].dy;
149 
150  // Calcul du Hessien
151  Warp->dWarpCompo(X1, X2, p, ptTemplateCompo[point].dW, dW);
152 
153  for (unsigned int it = 0; it < nbParam; it++)
154  tempt[it] = dW[0][it] * dIWx + dW[1][it] * dIWy;
155 
156  for (unsigned int it = 0; it < nbParam; it++)
157  for (unsigned int jt = 0; jt < nbParam; jt++)
158  HDir[it][jt] += tempt[it] * tempt[jt];
159 
160  for (unsigned int it = 0; it < nbParam; it++)
161  GDir[it] += er * tempt[it];
162  }
163  }
164  if (Nbpoint == 0) {
165  delete[] tempt;
166  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
167  }
168 
170 
171  try {
172  dp = (HLMDir).inverseByLU() * (GDir);
173  }
174  catch (const vpException &e) {
175  delete[] tempt;
176  throw(e);
177  }
178 
179  dp = gain * dp;
180  if (useBrent) {
181  alpha = 2.;
182  computeOptimalBrentGain(I, p, erreur / Nbpoint, dp, alpha);
183  dp = alpha * dp;
184  }
185 
186  p += dp;
187 
188  computeEvalRMS(p);
189 
190  if (iteration == 0) {
191  evolRMS_init = evolRMS;
192  }
193 
194  iteration++;
195 
196  evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
197  evolRMS_prec = evolRMS;
198 
199  } while ((iteration < iterationMax) && (evolRMS_delta > std::fabs(evolRMS_init) * evolRMS_eps));
200  delete[] tempt;
201 
202  nbIteration = iteration;
203 }
204 END_VISP_NAMESPACE
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:362
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:1143
error that can be emitted by ViSP classes.
Definition: vpException.h:60
@ badValue
Used to indicate that a value is not in the allowed range.
Definition: vpException.h:73
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
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:1873
void trackNoPyr(const vpImage< unsigned char > &I)
void initHessienDesired(const vpImage< unsigned char > &I)
VP_EXPLICIT vpTemplateTrackerSSDESM(vpTemplateTrackerWarp *warp)
void initCompInverse(const vpImage< unsigned char > &I)
virtual void dWarpCompo(const vpColVector &X1, const vpColVector &X2, const vpColVector &p, const double *dwdp0, vpMatrix &dM)=0
virtual bool isESMcompatible() const =0
virtual void getdW0(const int &v, const int &u, const double &dv, const double &du, double *dIdW)=0
virtual void getdWdp0(const int &v, const int &u, double *dIdW)=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
vpImage< double > BI
unsigned int templateSize
vpTemplateTrackerPointCompo * ptTemplateCompo
Error that can be emitted by the vpTracker class and its derivatives.
@ notEnoughPointError
Not enough point to track.