Visual Servoing Platform  version 3.6.1 under development (2024-12-07)
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 
44 BEGIN_VISP_NAMESPACE
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 {
72  initCompInverse(I);
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  moyIref = moyIref / Nbpoint;
145  moydIrefdp = moydIrefdp / Nbpoint;
146  moyd2Iref = moyd2Iref / Nbpoint;
147  moyIc = moyIc / Nbpoint;
148  Hdesire = 0;
149  double covarIref = 0, covarIc = 0;
150  double sIcIref = 0;
151  vpColVector sIcdIref(nbParam);
152  sIcdIref = 0;
153  vpMatrix sIcd2Iref(nbParam, nbParam);
154  sIcd2Iref = 0;
155  vpMatrix sdIrefdIref(nbParam, nbParam);
156  sdIrefdIref = 0;
157  for (unsigned int point = 0; point < templateSize; point++) {
158  i = ptTemplate[point].y;
159  j = ptTemplate[point].x;
160  X1[0] = j;
161  X1[1] = i;
162  X2[0] = j;
163  X2[1] = i;
164 
165  Warp->computeDenom(X1, p);
166 
167  j2 = X2[0];
168  i2 = X2[1];
169 
170  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
171  Iref = ptTemplate[point].val;
172 
173  if (!blur)
174  Ic = I.getValue(i2, j2);
175  else
176  Ic = BI.getValue(i2, j2);
177 
178  dIcx = dIx.getValue(i2, j2);
179  dIcy = dIy.getValue(i2, j2);
180 
181  Warp->dWarp(X1, X2, p, dW);
182 
183  for (unsigned int it = 0; it < nbParam; it++) {
184  tempt[it] = dW[0][it] * dIcx + dW[1][it] * dIcy;
185  }
186 
187  double prodIc = (Ic - moyIc);
188 
189  double d_Ixx = dIxx.getValue(i2, j2);
190  double d_Iyy = dIyy.getValue(i2, j2);
191  double d_Ixy = dIxy.getValue(i2, j2);
192 
193  for (unsigned int it = 0; it < nbParam; it++) {
194  for (unsigned int jt = 0; jt < nbParam; jt++) {
195  sIcd2Iref[it][jt] += prodIc * (dW[0][it] * (dW[0][jt] * d_Ixx + dW[1][jt] * d_Ixy) +
196  dW[1][it] * (dW[0][jt] * d_Ixy + dW[1][jt] * d_Iyy) - moyd2Iref[it][jt]);
197  sdIrefdIref[it][jt] +=
198  (ptTemplate[point].dW[it] - moydIrefdp[it]) * (ptTemplate[point].dW[jt] - moydIrefdp[jt]);
199  }
200  }
201 
202  for (unsigned int it = 0; it < nbParam; it++)
203  sIcdIref[it] += prodIc * (ptTemplate[point].dW[it] - moydIrefdp[it]);
204 
205  covarIref += (Iref - moyIref) * (Iref - moyIref);
206  covarIc += (Ic - moyIc) * (Ic - moyIc);
207  sIcIref += (Iref - moyIref) * (Ic - moyIc);
208  }
209  }
210  delete[] tempt;
211 
212  covarIref = sqrt(covarIref);
213  covarIc = sqrt(covarIc);
214 
215  denom = covarIref * covarIc;
216 
217  double NCC = sIcIref / denom;
218  vpColVector dcovarIref(nbParam);
219  dcovarIref = -sIcdIref / covarIref;
220 
221  vpColVector dNCC(nbParam);
222  dNCC = (sIcdIref / denom - NCC * dcovarIref / covarIref);
223  vpMatrix d2covarIref(nbParam, nbParam);
224  d2covarIref = -(sIcd2Iref - sdIrefdIref + dcovarIref * dcovarIref.t()) / covarIref;
225 #ifdef APPROX_NCC
226  Hdesire = sIcd2Iref / denom;
227 #else
228  Hdesire = (sIcd2Iref - sdIrefdIref + dcovarIref * dcovarIref.t()) / denom;
229 #endif
232 }
233 
235 {
236  if (blur)
238 
239  vpColVector dpinv(nbParam);
240  double Ic;
241  double Iref;
242  unsigned int iteration = 0;
243  int i, j;
244  double i2, j2;
245  initPosEvalRMS(p);
246 
247  double evolRMS_init = 0;
248  double evolRMS_prec = 0;
249  double evolRMS_delta;
250 
251  do {
252  unsigned int Nbpoint = 0;
253  G = 0;
254  Warp->computeCoeff(p);
255  double moyIref = 0;
256  double moyIc = 0;
257  for (unsigned int point = 0; point < templateSize; point++) {
258  i = ptTemplate[point].y;
259  j = ptTemplate[point].x;
260  X1[0] = j;
261  X1[1] = i;
262 
263  Warp->computeDenom(X1, p);
264  Warp->warpX(X1, X2, p);
265 
266  j2 = X2[0];
267  i2 = X2[1];
268  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
269  Iref = ptTemplate[point].val;
270 
271  if (!blur)
272  Ic = I.getValue(i2, j2);
273  else
274  Ic = BI.getValue(i2, j2);
275 
276  Nbpoint++;
277  moyIref += Iref;
278  moyIc += Ic;
279  }
280  }
281  if (Nbpoint > 0) {
282  moyIref = moyIref / Nbpoint;
283  moyIc = moyIc / Nbpoint;
284  double sIcIref = 0;
285  double covarIref = 0, covarIc = 0;
286  vpColVector sIcdIref(nbParam);
287  sIcdIref = 0;
288  vpColVector sIrefdIref(nbParam);
289  sIrefdIref = 0;
290 
291  for (unsigned int point = 0; point < templateSize; point++) {
292  i = ptTemplate[point].y;
293  j = ptTemplate[point].x;
294  X1[0] = j;
295  X1[1] = i;
296 
297  Warp->computeDenom(X1, p);
298  Warp->warpX(X1, X2, p);
299 
300  j2 = X2[0];
301  i2 = X2[1];
302  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
303  Iref = ptTemplate[point].val;
304 
305  if (!blur)
306  Ic = I.getValue(i2, j2);
307  else
308  Ic = BI.getValue(i2, j2);
309 
310  double prod = (Ic - moyIc);
311  for (unsigned int it = 0; it < nbParam; it++)
312  sIcdIref[it] += prod * (ptTemplate[point].dW[it] - moydIrefdp[it]);
313  for (unsigned int it = 0; it < nbParam; it++)
314  sIrefdIref[it] += (Iref - moyIref) * (ptTemplate[point].dW[it] - moydIrefdp[it]);
315 
316  covarIref += (Iref - moyIref) * (Iref - moyIref);
317  covarIc += (Ic - moyIc) * (Ic - moyIc);
318  sIcIref += (Iref - moyIref) * (Ic - moyIc);
319  }
320  }
321  covarIref = sqrt(covarIref);
322  covarIc = sqrt(covarIc);
323  double denom = covarIref * covarIc;
324 
325  if (std::fabs(denom) <= std::numeric_limits<double>::epsilon()) {
326  diverge = true;
327  }
328  else {
329  double NCC = sIcIref / denom;
330  vpColVector dcovarIref(nbParam);
331  dcovarIref = sIrefdIref / covarIref;
332  G = (sIcdIref / denom - NCC * dcovarIref / covarIref);
333 
334  try {
335  dp = -HLMdesireInverse * G;
336  }
337  catch (const vpException &e) {
338  throw(e);
339  }
340 
341  Warp->getParamInverse(dp, dpinv);
342  Warp->pRondp(p, dpinv, p);
343 
344  computeEvalRMS(p);
345  }
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 }
362 END_VISP_NAMESPACE
Implementation of column vector and the associated operations.
Definition: vpColVector.h:191
vpRowVector t() const
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:1143
error that can be emitted by ViSP classes.
Definition: vpException.h:60
static void getGradX(const vpImage< unsigned char > &I, vpImage< FilterType > &dIx, const vpImage< bool > *p_mask=nullptr)
static void getGradXGauss2D(const vpImage< ImageType > &I, vpImage< FilterType > &dIx, const FilterType *gaussianKernel, const FilterType *gaussianDerivativeKernel, unsigned int size, const vpImage< bool > *p_mask=nullptr)
static void filter(const vpImage< ImageType > &I, vpImage< FilterType > &If, const vpArray2D< FilterType > &M, bool convolve=false, const vpImage< bool > *p_mask=nullptr)
static void getGradYGauss2D(const vpImage< ImageType > &I, vpImage< FilterType > &dIy, const FilterType *gaussianKernel, const FilterType *gaussianDerivativeKernel, unsigned int size, const vpImage< bool > *p_mask=nullptr)
static void getGradY(const vpImage< unsigned char > &I, vpImage< FilterType > &dIy, 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:1869
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
VP_EXPLICIT vpTemplateTrackerZNCCInverseCompositional(vpTemplateTrackerWarp *warp)
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