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