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