Visual Servoing Platform version 3.5.0
vpTemplateTrackerSSDESM.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#include <visp3/core/vpImageFilter.h>
41#include <visp3/tt/vpTemplateTrackerSSDESM.h>
42
44 : vpTemplateTrackerSSD(warp), compoInitialised(false), HDir(), HInv(), HLMDir(), HLMInv(), GDir(), GInv()
45{
46 useCompositionnal = false;
47 useInverse = false;
48
49 if (!Warp->isESMcompatible()) {
51 "The selected warp function is not appropriate for the ESM algorithm..."));
52 }
53
60}
61
63
65{
67 int i, j;
68 // direct
69 for (unsigned int point = 0; point < templateSize; point++) {
70 i = ptTemplate[point].y;
71 j = ptTemplate[point].x;
72 ptTemplateCompo[point].dW = new double[2 * nbParam];
73 Warp->getdWdp0(i, j, ptTemplateCompo[point].dW);
74 }
75
76 // inverse
77 HInv = 0;
78 for (unsigned int point = 0; point < templateSize; point++) {
79 i = ptTemplate[point].y;
80 j = ptTemplate[point].x;
81
82 ptTemplate[point].dW = new double[nbParam];
83 Warp->getdW0(i, j, ptTemplate[point].dy, ptTemplate[point].dx, ptTemplate[point].dW);
84
85 for (unsigned int it = 0; it < nbParam; it++)
86 for (unsigned int jt = 0; jt < nbParam; jt++)
87 HInv[it][jt] += ptTemplate[point].dW[it] * ptTemplate[point].dW[jt];
88 }
90
91 compoInitialised = true;
92}
93
95{
96 if (blur)
100
101 double IW, dIWx, dIWy;
102 double Tij;
103 unsigned int iteration = 0;
104 int i, j;
105 double i2, j2;
106 double alpha = 2.;
107
109
110 double evolRMS_init = 0;
111 double evolRMS_prec = 0;
112 double evolRMS_delta;
113 double *tempt = new double[nbParam];
114
115 do {
116 unsigned int Nbpoint = 0;
117 double erreur = 0;
118 dp = 0;
119 HDir = 0;
120 GDir = 0;
121 GInv = 0;
122 Warp->computeCoeff(p);
123 for (unsigned int point = 0; point < templateSize; point++) {
124 i = ptTemplate[point].y;
125 j = ptTemplate[point].x;
126 X1[0] = j;
127 X1[1] = i;
128
129 Warp->computeDenom(X1, p);
130 Warp->warpX(X1, X2, p);
131
132 j2 = X2[0];
133 i2 = X2[1];
134 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
135 // INVERSE
136 Tij = ptTemplate[point].val;
137 if (!blur)
138 IW = I.getValue(i2, j2);
139 else
140 IW = BI.getValue(i2, j2);
141 Nbpoint++;
142 double er = (Tij - IW);
143 for (unsigned int it = 0; it < nbParam; it++)
144 GInv[it] += er * ptTemplate[point].dW[it];
145
146 erreur += er * er;
147
148 dIWx = dIx.getValue(i2, j2) + ptTemplate[point].dx;
149 dIWy = dIy.getValue(i2, j2) + ptTemplate[point].dy;
150
151 // Calcul du Hessien
152 Warp->dWarpCompo(X1, X2, p, ptTemplateCompo[point].dW, dW);
153
154 for (unsigned int it = 0; it < nbParam; it++)
155 tempt[it] = dW[0][it] * dIWx + dW[1][it] * dIWy;
156
157 for (unsigned int it = 0; it < nbParam; it++)
158 for (unsigned int jt = 0; jt < nbParam; jt++)
159 HDir[it][jt] += tempt[it] * tempt[jt];
160
161 for (unsigned int it = 0; it < nbParam; it++)
162 GDir[it] += er * tempt[it];
163 }
164 }
165 if (Nbpoint == 0) {
166 delete[] tempt;
167 throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
168 }
169
171
172 try {
173 dp = (HLMDir).inverseByLU() * (GDir);
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
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}
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:304
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:310
error that can be emited by ViSP classes.
Definition: vpException.h:72
@ badValue
Used to indicate that a value is not in the allowed range.
Definition: vpException.h:97
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
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:6683
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 emited by the vpTracker class and its derivates.