Visual Servoing Platform  version 3.3.1 under development (2020-10-22)
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 {
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 
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;
247  initPosEvalRMS(p);
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 
345  computeEvalRMS(p);
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 }
virtual void warpX(const int &v1, const int &u1, double &v2, double &u2, const vpColVector &p)=0
virtual void getParamInverse(const vpColVector &p, vpColVector &p_inv) const =0
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:156
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1346
static void getGradX(const vpImage< unsigned char > &I, vpImage< double > &dIx)
vpTemplateTrackerPoint * ptTemplate
static void getGradY(const vpImage< unsigned char > &I, vpImage< double > &dIy)
vpMatrix inverseByLU() const
error that can be emited by ViSP classes.
Definition: vpException.h:71
void initPosEvalRMS(const vpColVector &p)
vpRowVector t() const
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
void computeEvalRMS(const vpColVector &p)
static void getGradXGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIx, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
vpImage< double > BI
unsigned int templateSize
unsigned int iterationMax
vpImage< double > dIx
virtual void pRondp(const vpColVector &p1, const vpColVector &p2, vpColVector &p12) const =0
vpImage< double > dIy
virtual void getdW0(const int &v, const int &u, const double &dv, const double &du, double *dIdW)=0
unsigned int nbIteration
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, bool convolve=false)
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:310
unsigned int getHeight() const
Definition: vpImage.h:188
Implementation of column vector and the associated operations.
Definition: vpColVector.h:130
virtual void dWarp(const vpColVector &X1, const vpColVector &X2, const vpColVector &p, vpMatrix &dM)=0
vpTemplateTrackerWarp * Warp
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:6779
unsigned int getWidth() const
Definition: vpImage.h:246