Visual Servoing Platform  version 3.6.1 under development (2024-11-07)
vpTemplateTrackerSSDInverseCompositional.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/vpImageTools.h>
40 #include <visp3/tt/vpTemplateTrackerSSDInverseCompositional.h>
41 
42 BEGIN_VISP_NAMESPACE
44  : vpTemplateTrackerSSD(warp), compoInitialised(false), HInv(), HCompInverse(), useTemplateSelect(false)
45 {
46  useInverse = true;
49 }
50 
52 {
53  H = 0;
54 
55  for (unsigned int point = 0; point < templateSize; point++) {
56  if ((!useTemplateSelect) || (ptTemplateSelect[point])) {
57  int i = ptTemplate[point].y;
58  int j = ptTemplate[point].x;
59  ptTemplate[point].dW = new double[nbParam];
60 
61  Warp->getdW0(i, j, ptTemplate[point].dy, ptTemplate[point].dx, ptTemplate[point].dW);
62 
63  for (unsigned int it = 0; it < nbParam; it++)
64  for (unsigned int jt = 0; jt < nbParam; jt++)
65  H[it][jt] += ptTemplate[point].dW[it] * ptTemplate[point].dW[jt];
66  }
67  }
68 
69  vpMatrix HLMtemp;
70  vpMatrix::computeHLM(H, lambdaDep, HLMtemp);
71 
72  HCompInverse = HLMtemp.inverseByLU();
73 
74  for (unsigned int point = 0; point < templateSize; point++) {
75  if ((!useTemplateSelect) || (ptTemplateSelect[point])) {
76  ptTemplate[point].HiG = new double[nbParam];
77 
78  for (unsigned int i = 0; i < HCompInverse.getRows(); i++) {
79  ptTemplate[point].HiG[i] = 0;
80  for (unsigned int j = 0; j < HCompInverse.getCols(); j++) {
81  ptTemplate[point].HiG[i] -= HCompInverse[i][j] * ptTemplate[point].dW[j];
82  }
83  }
84  }
85  }
86  compoInitialised = true;
87 }
88 
90 {
91  initCompInverse(I);
92 }
93 
95 {
96  if (blur)
98 
99  vpColVector dpinv(nbParam);
100  double IW;
101  double Tij;
102  unsigned int iteration = 0;
103  int i, j;
104  double i2, j2;
105  double alpha = 2.;
106  initPosEvalRMS(p);
107 
109 
110  double evolRMS_init = 0;
111  double evolRMS_prec = 0;
112  double evolRMS_delta;
113 
114  do {
115  unsigned int Nbpoint = 0;
116  double erreur = 0;
117  dp = 0;
118  Warp->computeCoeff(p);
119  for (unsigned int point = 0; point < templateSize; point++) {
120  if ((!useTemplateSelect) || (ptTemplateSelect[point])) {
121  pt = &ptTemplate[point];
122  i = pt->y;
123  j = pt->x;
124  X1[0] = j;
125  X1[1] = i;
126  Warp->computeDenom(X1, p);
127  Warp->warpX(X1, X2, p);
128  j2 = X2[0];
129  i2 = X2[1];
130 
131  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
132  Tij = pt->val;
133  if (!blur)
134  IW = I.getValue(i2, j2);
135  else
136  IW = BI.getValue(i2, j2);
137  Nbpoint++;
138  double er = (Tij - IW);
139  for (unsigned int it = 0; it < nbParam; it++)
140  dp[it] += er * pt->HiG[it];
141 
142  erreur += er * er;
143  }
144  }
145  }
146  if (Nbpoint == 0) {
147  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
148  }
149  dp = gain * dp;
150  if (useBrent) {
151  alpha = 2.;
152  computeOptimalBrentGain(I, p, erreur / Nbpoint, dp, alpha);
153  dp = alpha * dp;
154  }
155  Warp->getParamInverse(dp, dpinv);
156  Warp->pRondp(p, dpinv, p);
157 
158  computeEvalRMS(p);
159 
160  if (iteration == 0) {
161  evolRMS_init = evolRMS;
162  }
163  iteration++;
164 
165  evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
166  evolRMS_prec = evolRMS;
167 
168  } while ((iteration < iterationMax) && (evolRMS_delta > std::fabs(evolRMS_init) * evolRMS_eps));
169  nbIteration = iteration;
170 }
171 END_VISP_NAMESPACE
unsigned int getCols() const
Definition: vpArray2D.h:337
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:362
unsigned int getRows() const
Definition: vpArray2D.h:347
Implementation of column vector and the associated operations.
Definition: vpColVector.h:191
static void filter(const vpImage< ImageType > &I, vpImage< FilterType > &If, const vpArray2D< FilterType > &M, bool convolve=false, 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
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:169
vpMatrix inverseByLU() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:1873
VP_EXPLICIT vpTemplateTrackerSSDInverseCompositional(vpTemplateTrackerWarp *warp)
virtual void getParamInverse(const vpColVector &p, vpColVector &p_inv) const =0
virtual void getdW0(const int &v, const int &u, const double &dv, const double &du, double *dIdW)=0
virtual void warpX(const int &v1, const int &u1, double &v2, double &u2, const vpColVector &p)=0
virtual void pRondp(const vpColVector &p1, const vpColVector &p2, vpColVector &p12) const =0
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
Error that can be emitted by the vpTracker class and its derivatives.
@ notEnoughPointError
Not enough point to track.