Visual Servoing Platform version 3.5.0
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
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 } catch (const vpException &e) {
129 delete[] tempt;
130 throw(e);
131 }
132
133 switch (minimizationMethod) {
135 vpColVector p_test_LMA(nbParam);
136 p_test_LMA = p + dp;
137 erreur = -getCost(I, p);
138 double erreur_LMA = -getCost(I, p_test_LMA);
139 if (erreur_LMA < erreur) {
140 p = p_test_LMA;
141 lambda = (lambda / 10. < 1e-6) ? lambda / 10. : 1e-6;
142 } else {
143 lambda = (lambda * 10. < 1e6) ? 1e6 : lambda * 10.;
144 }
145 } break;
147 dp = gain * 0.000001 * G;
148 if (useBrent) {
149 alpha = 2.;
150 computeOptimalBrentGain(I, p, erreur, dp, alpha);
151 dp = alpha * dp;
152 }
153 p += dp;
154 break;
155 }
156
158 if (iterationGlobale != 0) {
159 vpColVector s_quasi = p - p_prec;
160 vpColVector y_quasi = G - G_prec;
161 double s_scal_y = s_quasi.t() * y_quasi;
162 if (std::fabs(s_scal_y) > std::numeric_limits<double>::epsilon()) // DFP
163 KQuasiNewton = KQuasiNewton + 0.001 * (s_quasi * s_quasi.t() / s_scal_y -
164 KQuasiNewton * y_quasi * y_quasi.t() * KQuasiNewton /
165 (y_quasi.t() * KQuasiNewton * y_quasi));
166 }
167 dp = -KQuasiNewton * G;
168 p_prec = p;
169 G_prec = G;
170 p -= 1.01 * dp;
171 } break;
172
174 default: {
175 if (useBrent) {
176 alpha = 2.;
177 computeOptimalBrentGain(I, p, erreur, dp, alpha);
178 dp = alpha * dp;
179 }
180
181 p += dp;
182 break;
183 }
184 }
185
187
188 if (iteration == 0) {
189 evolRMS_init = evolRMS;
190 }
191
192 iteration++;
194
195 evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
196 evolRMS_prec = evolRMS;
197
198 } while ( (iteration < iterationMax) && (evolRMS_delta > std::fabs(evolRMS_init)*evolRMS_eps) );
199 delete[] tempt;
200
201 nbIteration = iteration;
202}
Implementation of column vector and the associated operations.
Definition: vpColVector.h:131
vpRowVector t() const
error that can be emited by ViSP classes.
Definition: vpException.h:72
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, bool convolve=false)
static void getGradXGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIx, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
unsigned int getWidth() const
Definition: vpImage.h:246
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1346
unsigned int getHeight() const
Definition: vpImage.h:188
vpMatrix inverseByLU() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:6683
void trackNoPyr(const vpImage< unsigned char > &I)
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 emited by the vpTracker class and its derivates.