Visual Servoing Platform  version 3.6.1 under development (2024-03-28)
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 
45  : vpTemplateTrackerZNCC(warp), compoInitialised(false), moydIrefdp()
46 {
47  useInverse = true;
48 }
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 {
71  initCompInverse(I);
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;
244  initPosEvalRMS(p);
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 
341  computeEvalRMS(p);
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.
Definition: vpColVector.h:163
vpRowVector t() const
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:1056
error that can be emitted by ViSP classes.
Definition: vpException.h:59
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:245
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1764
unsigned int getHeight() const
Definition: vpImage.h:184
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:146
vpMatrix inverseByLU() const
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:6407
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