Visual Servoing Platform version 3.5.0
vpTemplateTrackerZNCCInverseCompositional.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 <limits> // numeric_limits
41
42#include <visp3/core/vpImageFilter.h>
43#include <visp3/tt/vpTemplateTrackerZNCCInverseCompositional.h>
44
46 : vpTemplateTrackerZNCC(warp), compoInitialised(false), moydIrefdp()
47{
48 useInverse = true;
49}
50
52{
55
56 for (unsigned int point = 0; point < templateSize; point++) {
57 int i = ptTemplate[point].y;
58 int j = ptTemplate[point].x;
59
60 ptTemplate[point].dW = new double[nbParam];
61
62 double dx = ptTemplate[point].dx;
63 double dy = ptTemplate[point].dy;
64
65 Warp->getdW0(i, j, dy, dx, ptTemplate[point].dW);
66 }
67 compoInitialised = true;
68}
69
71{
73
74 if (blur)
78
79 vpImage<double> dIxx, dIxy, dIyx, dIyy;
82
85
86 Warp->computeCoeff(p);
87 double Ic, dIcx = 0., dIcy = 0.;
88 double Iref;
89 int i, j;
90 double i2, j2;
91 int Nbpoint = 0;
92
93 double moyIref = 0;
94 double moyIc = 0;
95 double denom = 0;
96 double *tempt = new double[nbParam];
97
99 moydIrefdp = 0;
100 vpMatrix moyd2Iref(nbParam, nbParam);
101 moyd2Iref = 0;
102
103 for (unsigned int point = 0; point < templateSize; point++) {
104 i = ptTemplate[point].y;
105 j = ptTemplate[point].x;
106 X1[0] = j;
107 X1[1] = i;
108 X2[0] = j;
109 X2[1] = i;
110
111 j2 = X2[0];
112 i2 = X2[1];
113
114 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
115 Iref = ptTemplate[point].val;
116
117 if (!blur)
118 Ic = I.getValue(i2, j2);
119 else
120 Ic = BI.getValue(i2, j2);
121
122 Nbpoint++;
123 moyIref += Iref;
124 moyIc += Ic;
125
126 for (unsigned int it = 0; it < nbParam; it++)
127 moydIrefdp[it] += ptTemplate[point].dW[it];
128
129 Warp->dWarp(X1, X2, p, dW);
130 for (unsigned int it = 0; it < nbParam; it++)
131 tempt[it] = dW[0][it] * dIcx + dW[1][it] * dIcy;
132 double d_Ixx = dIxx.getValue(i2, j2);
133 double d_Iyy = dIyy.getValue(i2, j2);
134 double d_Ixy = dIxy.getValue(i2, j2);
135
136 for (unsigned int it = 0; it < nbParam; it++)
137 for (unsigned int jt = 0; jt < nbParam; jt++) {
138 moyd2Iref[it][jt] += (dW[0][it] * (dW[0][jt] * d_Ixx + dW[1][jt] * d_Ixy) +
139 dW[1][it] * (dW[0][jt] * d_Ixy + dW[1][jt] * d_Iyy));
140 }
141
142 }
143 }
144
145 moyIref = moyIref / Nbpoint;
146 moydIrefdp = moydIrefdp / Nbpoint;
147 moyd2Iref = moyd2Iref / Nbpoint;
148 moyIc = moyIc / Nbpoint;
149 Hdesire = 0;
150 double covarIref = 0, covarIc = 0;
151 double sIcIref = 0;
152 vpColVector sIcdIref(nbParam);
153 sIcdIref = 0;
154 vpMatrix sIcd2Iref(nbParam, nbParam);
155 sIcd2Iref = 0;
156 vpMatrix sdIrefdIref(nbParam, nbParam);
157 sdIrefdIref = 0;
158 for (unsigned int point = 0; point < templateSize; point++) {
159 i = ptTemplate[point].y;
160 j = ptTemplate[point].x;
161 X1[0] = j;
162 X1[1] = i;
163 X2[0] = j;
164 X2[1] = i;
165
166 Warp->computeDenom(X1, p);
167
168 j2 = X2[0];
169 i2 = X2[1];
170
171 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
172 Iref = ptTemplate[point].val;
173
174 if (!blur)
175 Ic = I.getValue(i2, j2);
176 else
177 Ic = BI.getValue(i2, j2);
178
179 dIcx = dIx.getValue(i2, j2);
180 dIcy = dIy.getValue(i2, j2);
181
182 Warp->dWarp(X1, X2, p, dW);
183
184 for (unsigned int it = 0; it < nbParam; it++) {
185 tempt[it] = dW[0][it] * dIcx + dW[1][it] * dIcy;
186 }
187
188 double prodIc = (Ic - moyIc);
189
190 double d_Ixx = dIxx.getValue(i2, j2);
191 double d_Iyy = dIyy.getValue(i2, j2);
192 double d_Ixy = dIxy.getValue(i2, j2);
193
194 for (unsigned int it = 0; it < nbParam; it++) {
195 for (unsigned int jt = 0; jt < nbParam; jt++) {
196 sIcd2Iref[it][jt] += prodIc * (dW[0][it] * (dW[0][jt] * d_Ixx + dW[1][jt] * d_Ixy) +
197 dW[1][it] * (dW[0][jt] * d_Ixy + dW[1][jt] * d_Iyy) - moyd2Iref[it][jt]);
198 sdIrefdIref[it][jt] +=
199 (ptTemplate[point].dW[it] - moydIrefdp[it]) * (ptTemplate[point].dW[jt] - moydIrefdp[jt]);
200 }
201 }
202
203
204 for (unsigned int it = 0; it < nbParam; it++)
205 sIcdIref[it] += prodIc * (ptTemplate[point].dW[it] - moydIrefdp[it]);
206
207 covarIref += (Iref - moyIref) * (Iref - moyIref);
208 covarIc += (Ic - moyIc) * (Ic - moyIc);
209 sIcIref += (Iref - moyIref) * (Ic - moyIc);
210 }
211 }
212 delete[] tempt;
213
214 covarIref = sqrt(covarIref);
215 covarIc = sqrt(covarIc);
216
217 denom = covarIref * covarIc;
218
219 double NCC = sIcIref / denom;
220 vpColVector dcovarIref(nbParam);
221 dcovarIref = -sIcdIref / covarIref;
222
223 vpColVector dNCC(nbParam);
224 dNCC = (sIcdIref / denom - NCC * dcovarIref / covarIref);
225 vpMatrix d2covarIref(nbParam, nbParam);
226 d2covarIref = -(sIcd2Iref - sdIrefdIref + dcovarIref * dcovarIref.t()) / covarIref;
227#ifdef APPROX_NCC
228 Hdesire = sIcd2Iref / denom;
229#else
230 Hdesire = (sIcd2Iref - sdIrefdIref + dcovarIref * dcovarIref.t()) / denom;
231#endif
234}
235
237{
238 if (blur)
240
241 vpColVector dpinv(nbParam);
242 double Ic;
243 double Iref;
244 unsigned int iteration = 0;
245 int i, j;
246 double i2, j2;
248
249 double evolRMS_init = 0;
250 double evolRMS_prec = 0;
251 double evolRMS_delta;
252
253 do {
254 unsigned int Nbpoint = 0;
255 G = 0;
256 Warp->computeCoeff(p);
257 double moyIref = 0;
258 double moyIc = 0;
259 for (unsigned int point = 0; point < templateSize; point++) {
260 i = ptTemplate[point].y;
261 j = ptTemplate[point].x;
262 X1[0] = j;
263 X1[1] = i;
264
265 Warp->computeDenom(X1, p);
266 Warp->warpX(X1, X2, p);
267
268 j2 = X2[0];
269 i2 = X2[1];
270 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
271 Iref = ptTemplate[point].val;
272
273 if (!blur)
274 Ic = I.getValue(i2, j2);
275 else
276 Ic = BI.getValue(i2, j2);
277
278 Nbpoint++;
279 moyIref += Iref;
280 moyIc += Ic;
281 }
282 }
283 if (Nbpoint > 0) {
284 moyIref = moyIref / Nbpoint;
285 moyIc = moyIc / Nbpoint;
286 double sIcIref = 0;
287 double covarIref = 0, covarIc = 0;
288 vpColVector sIcdIref(nbParam);
289 sIcdIref = 0;
290 vpColVector sIrefdIref(nbParam);
291 sIrefdIref = 0;
292
293 for (unsigned int point = 0; point < templateSize; point++) {
294 i = ptTemplate[point].y;
295 j = ptTemplate[point].x;
296 X1[0] = j;
297 X1[1] = i;
298
299 Warp->computeDenom(X1, p);
300 Warp->warpX(X1, X2, p);
301
302 j2 = X2[0];
303 i2 = X2[1];
304 if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
305 Iref = ptTemplate[point].val;
306
307 if (!blur)
308 Ic = I.getValue(i2, j2);
309 else
310 Ic = BI.getValue(i2, j2);
311
312 double prod = (Ic - moyIc);
313 for (unsigned int it = 0; it < nbParam; it++)
314 sIcdIref[it] += prod * (ptTemplate[point].dW[it] - moydIrefdp[it]);
315 for (unsigned int it = 0; it < nbParam; it++)
316 sIrefdIref[it] += (Iref - moyIref) * (ptTemplate[point].dW[it] - moydIrefdp[it]);
317
318 covarIref += (Iref - moyIref) * (Iref - moyIref);
319 covarIc += (Ic - moyIc) * (Ic - moyIc);
320 sIcIref += (Iref - moyIref) * (Ic - moyIc);
321 }
322 }
323 covarIref = sqrt(covarIref);
324 covarIc = sqrt(covarIc);
325 double denom = covarIref * covarIc;
326
327 if (std::fabs(denom) <= std::numeric_limits<double>::epsilon()) {
328 diverge = true;
329 } else {
330 double NCC = sIcIref / denom;
331 vpColVector dcovarIref(nbParam);
332 dcovarIref = sIrefdIref / covarIref;
333 G = (sIcdIref / denom - NCC * dcovarIref / covarIref);
334
335 try {
336 dp = - HLMdesireInverse * G;
337 }
338 catch (const vpException &e) {
339 throw(e);
340 }
341
342 Warp->getParamInverse(dp, dpinv);
343 Warp->pRondp(p, dpinv, p);
344
346 }
347 } else
348 diverge = true;
349
350 if (iteration == 0) {
351 evolRMS_init = evolRMS;
352 }
353 iteration++;
354
355 evolRMS_delta = std::fabs(evolRMS - evolRMS_prec);
356 evolRMS_prec = evolRMS;
357
358 } while ( (!diverge && (evolRMS_delta > std::fabs(evolRMS_init)*evolRMS_eps) && (iteration < iterationMax)) );
359
360 nbIteration = iteration;
361}
Implementation of column vector and the associated operations.
Definition: vpColVector.h:131
vpRowVector t() const
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:310
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 getGradX(const vpImage< unsigned char > &I, vpImage< double > &dIx)
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
static void getGradY(const vpImage< unsigned char > &I, vpImage< double > &dIy)
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
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:154
vpMatrix inverseByLU() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:6683
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 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
virtual void pRondp(const vpColVector &p1, const vpColVector &p2, vpColVector &p12) const =0
vpImage< double > dIx
vpImage< double > dIy
void computeEvalRMS(const vpColVector &p)
unsigned int iterationMax
void initPosEvalRMS(const vpColVector &p)
unsigned int nbIteration
vpTemplateTrackerPoint * ptTemplate
vpTemplateTrackerWarp * Warp
vpImage< double > BI
unsigned int templateSize