Visual Servoing Platform  version 3.2.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), evolRMS(0), x_pos(), y_pos(), threshold_RMS(1e-8),
47  moydIrefdp()
48 {
49  useInverse = true;
50 }
51 
53 {
54  // std::cout<<"Initialise precomputed value of Compositionnal
55  // Inverse"<<std::endl;
58 
59  for (unsigned int point = 0; point < templateSize; point++) {
60  int i = ptTemplate[point].y;
61  int j = ptTemplate[point].x;
62 
63  X1[0] = j;
64  X1[1] = i;
65  Warp->computeDenom(X1, p);
66  ptTemplate[point].dW = new double[nbParam];
67 
68  double dx = ptTemplate[point].dx;
69  double dy = ptTemplate[point].dy;
70  // std::cout<<ptTemplate[point].dx<<","<<ptTemplate[point].dy<<std::endl;
71 
72  Warp->getdW0(i, j, dy, dx, ptTemplate[point].dW);
73  }
74  // vpTRACE("fin Comp Inverse");
75  compoInitialised = true;
76 }
77 
79 {
80  initCompInverse(I);
81 
82  if (blur)
86 
87  vpImage<double> dIxx, dIxy, dIyx, dIyy;
90 
93 
94  Warp->computeCoeff(p);
95  double Ic, dIcx = 0., dIcy = 0.;
96  double Iref;
97  int i, j;
98  double i2, j2;
99  int Nbpoint = 0;
100 
101  double moyIref = 0;
102  double moyIc = 0;
103  double denom = 0;
105  moydIrefdp = 0;
106  vpMatrix moyd2Iref(nbParam, nbParam);
107  moyd2Iref = 0;
108 
109  for (unsigned int point = 0; point < templateSize; point++) {
110  i = ptTemplate[point].y;
111  j = ptTemplate[point].x;
112  X1[0] = j;
113  X1[1] = i;
114  X2[0] = j;
115  X2[1] = i;
116 
117  Warp->computeDenom(X1, p);
118 
119  j2 = X2[0];
120  i2 = X2[1];
121 
122  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
123  Iref = ptTemplate[point].val;
124 
125  if (!blur)
126  Ic = I.getValue(i2, j2);
127  else
128  Ic = BI.getValue(i2, j2);
129 
130  Nbpoint++;
131  moyIref += Iref;
132  moyIc += Ic;
133 
134  for (unsigned int it = 0; it < nbParam; it++)
135  moydIrefdp[it] += ptTemplate[point].dW[it];
136 
137  Warp->dWarp(X1, X2, p, dW);
138  double *tempt = new double[nbParam];
139  for (unsigned int it = 0; it < nbParam; it++)
140  tempt[it] = dW[0][it] * dIcx + dW[1][it] * dIcy;
141  double d_Ixx = dIxx.getValue(i2, j2);
142  double d_Iyy = dIyy.getValue(i2, j2);
143  double d_Ixy = dIxy.getValue(i2, j2);
144 
145  for (unsigned int it = 0; it < nbParam; it++)
146  for (unsigned int jt = 0; jt < nbParam; jt++) {
147  moyd2Iref[it][jt] += (dW[0][it] * (dW[0][jt] * d_Ixx + dW[1][jt] * d_Ixy) +
148  dW[1][it] * (dW[0][jt] * d_Ixy + dW[1][jt] * d_Iyy));
149  }
150 
151  delete[] tempt;
152  }
153  }
154 
155  moyIref = moyIref / Nbpoint;
156  moydIrefdp = moydIrefdp / Nbpoint;
157  moyd2Iref = moyd2Iref / Nbpoint;
158  moyIc = moyIc / Nbpoint;
159  Hdesire = 0;
160  double covarIref = 0, covarIc = 0;
161  double sIcIref = 0;
162  vpColVector sIcdIref(nbParam);
163  sIcdIref = 0;
164  vpMatrix sIcd2Iref(nbParam, nbParam);
165  sIcd2Iref = 0;
166  vpMatrix sdIrefdIref(nbParam, nbParam);
167  sdIrefdIref = 0;
168  for (unsigned int point = 0; point < templateSize; point++) {
169  i = ptTemplate[point].y;
170  j = ptTemplate[point].x;
171  X1[0] = j;
172  X1[1] = i;
173  X2[0] = j;
174  X2[1] = i;
175 
176  Warp->computeDenom(X1, p);
177 
178  j2 = X2[0];
179  i2 = X2[1];
180 
181  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
182  Iref = ptTemplate[point].val;
183 
184  if (!blur)
185  Ic = I.getValue(i2, j2);
186  else
187  Ic = BI.getValue(i2, j2);
188 
189  dIcx = dIx.getValue(i2, j2);
190  dIcy = dIy.getValue(i2, j2);
191 
192  Warp->dWarp(X1, X2, p, dW);
193 
194  double *tempt = new double[nbParam];
195  for (unsigned int it = 0; it < nbParam; it++)
196  tempt[it] = dW[0][it] * dIcx + dW[1][it] * dIcy;
197 
198  double prodIc = (Ic - moyIc);
199 
200  double d_Ixx = dIxx.getValue(i2, j2);
201  double d_Iyy = dIyy.getValue(i2, j2);
202  double d_Ixy = dIxy.getValue(i2, j2);
203 
204  for (unsigned int it = 0; it < nbParam; it++)
205  for (unsigned int jt = 0; jt < nbParam; jt++) {
206  sIcd2Iref[it][jt] += prodIc * (dW[0][it] * (dW[0][jt] * d_Ixx + dW[1][jt] * d_Ixy) +
207  dW[1][it] * (dW[0][jt] * d_Ixy + dW[1][jt] * d_Iyy) - moyd2Iref[it][jt]);
208  sdIrefdIref[it][jt] +=
209  (ptTemplate[point].dW[it] - moydIrefdp[it]) * (ptTemplate[point].dW[jt] - moydIrefdp[jt]);
210  }
211 
212  delete[] tempt;
213 
214  for (unsigned int it = 0; it < nbParam; it++)
215  sIcdIref[it] += prodIc * (ptTemplate[point].dW[it] - moydIrefdp[it]);
216 
217  covarIref += (Iref - moyIref) * (Iref - moyIref);
218  covarIc += (Ic - moyIc) * (Ic - moyIc);
219  sIcIref += (Iref - moyIref) * (Ic - moyIc);
220  }
221  }
222  covarIref = sqrt(covarIref);
223  covarIc = sqrt(covarIc);
224 
225  denom = covarIref * covarIc;
226 
227  double NCC = sIcIref / denom;
228  // std::cout<<"NCC = "<<NCC<<std::endl;
229  vpColVector dcovarIref(nbParam);
230  dcovarIref = -sIcdIref / covarIref;
231 
232  vpColVector dNCC(nbParam);
233  dNCC = (sIcdIref / denom - NCC * dcovarIref / covarIref);
234  vpMatrix d2covarIref(nbParam, nbParam);
235  d2covarIref = -(sIcd2Iref - sdIrefdIref + dcovarIref * dcovarIref.t()) / covarIref;
236 #ifdef APPROX_NCC
237  Hdesire = sIcd2Iref / denom;
238 #else
239  Hdesire = (sIcd2Iref - sdIrefdIref + dcovarIref * dcovarIref.t()) / denom;
240 #endif
243  // std::cout<<"Hdesire = "<<Hdesire<<std::endl;
244 }
245 
247 {
248  if (blur)
250 
251  // double erreur=0;
252  vpColVector dpinv(nbParam);
253  double Ic;
254  double Iref;
255  unsigned int iteration = 0;
256  int i, j;
257  double i2, j2;
258  initPosEvalRMS(p);
259  do {
260  unsigned int Nbpoint = 0;
261  // erreur=0;
262  G = 0;
263  Warp->computeCoeff(p);
264  double moyIref = 0;
265  double moyIc = 0;
266  for (unsigned int point = 0; point < templateSize; point++) {
267  i = ptTemplate[point].y;
268  j = ptTemplate[point].x;
269  X1[0] = j;
270  X1[1] = i;
271 
272  Warp->computeDenom(X1, p);
273  Warp->warpX(X1, X2, p);
274 
275  j2 = X2[0];
276  i2 = X2[1];
277  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
278  Iref = ptTemplate[point].val;
279 
280  if (!blur)
281  Ic = I.getValue(i2, j2);
282  else
283  Ic = BI.getValue(i2, j2);
284 
285  Nbpoint++;
286  moyIref += Iref;
287  moyIc += Ic;
288  }
289  }
290  if (Nbpoint > 0) {
291  moyIref = moyIref / Nbpoint;
292  moyIc = moyIc / Nbpoint;
293  double sIcIref = 0;
294  double covarIref = 0, covarIc = 0;
295  vpColVector sIcdIref(nbParam);
296  sIcdIref = 0;
297  vpColVector sIrefdIref(nbParam);
298  sIrefdIref = 0;
299 
300  for (unsigned int point = 0; point < templateSize; point++) {
301  i = ptTemplate[point].y;
302  j = ptTemplate[point].x;
303  X1[0] = j;
304  X1[1] = i;
305 
306  Warp->computeDenom(X1, p);
307  Warp->warpX(X1, X2, p);
308 
309  j2 = X2[0];
310  i2 = X2[1];
311  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
312  Iref = ptTemplate[point].val;
313 
314  if (!blur)
315  Ic = I.getValue(i2, j2);
316  else
317  Ic = BI.getValue(i2, j2);
318 
319  double prod = (Ic - moyIc);
320  for (unsigned int it = 0; it < nbParam; it++)
321  sIcdIref[it] += prod * (ptTemplate[point].dW[it] - moydIrefdp[it]);
322  for (unsigned int it = 0; it < nbParam; it++)
323  sIrefdIref[it] += (Iref - moyIref) * (ptTemplate[point].dW[it] - moydIrefdp[it]);
324 
325  // double er=(Iref-Ic);
326  // erreur+=(er*er);
327  // denom+=(Iref-moyIref)*(Iref-moyIref)*(Ic-moyIc)*(Ic-moyIc);
328  covarIref += (Iref - moyIref) * (Iref - moyIref);
329  covarIc += (Ic - moyIc) * (Ic - moyIc);
330  sIcIref += (Iref - moyIref) * (Ic - moyIc);
331  }
332  }
333  covarIref = sqrt(covarIref);
334  covarIc = sqrt(covarIc);
335  double denom = covarIref * covarIc;
336 
337  // if(denom==0.0)
338  if (std::fabs(denom) <= std::numeric_limits<double>::epsilon()) {
339  diverge = true;
340  } else {
341  double NCC = sIcIref / denom;
342  vpColVector dcovarIref(nbParam);
343  dcovarIref = sIrefdIref / covarIref;
344  G = 1. * (sIcdIref / denom - NCC * dcovarIref / covarIref);
345 
346  try {
347  dp = -1. * HLMdesireInverse * G;
348  } catch (...) {
349  std::cout << "probleme inversion" << std::endl;
350  break;
351  }
352 
353  Warp->getParamInverse(dp, dpinv);
354  Warp->pRondp(p, dpinv, p);
355 
356  computeEvalRMS(p);
357  }
358  } else
359  diverge = true;
360 
361  iteration++;
362  } while ((!diverge && (evolRMS > threshold_RMS) && (iteration < iterationMax)));
363 
364  // std::cout<<"erreur "<<erreur<<std::endl;
365  nbIteration = iteration;
366 
368 }
369 
371 {
372  unsigned int nb_corners = zoneTracked->getNbTriangle() * 3;
373  x_pos.resize(nb_corners);
374  y_pos.resize(nb_corners);
375 
376  Warp->computeCoeff(p);
377  vpTemplateTrackerTriangle triangle;
378 
379  for (unsigned int i = 0; i < zoneTracked->getNbTriangle(); i++) {
380  zoneTracked->getTriangle(i, triangle);
381  for (unsigned int j = 0; j < 3; j++) {
382  triangle.getCorner(j, X1[0], X1[1]);
383 
384  Warp->computeDenom(X1, p_);
385  Warp->warpX(X1, X2, p_);
386  x_pos[i * 3 + j] = X2[0];
387  y_pos[i * 3 + j] = X2[1];
388  }
389  }
390 }
391 
393 {
394  unsigned int nb_corners = zoneTracked->getNbTriangle() * 3;
395 
396  Warp->computeCoeff(p_);
397  evolRMS = 0;
398  vpTemplateTrackerTriangle triangle;
399 
400  for (unsigned int i = 0; i < zoneTracked->getNbTriangle(); i++) {
401  zoneTracked->getTriangle(i, triangle);
402  for (unsigned int j = 0; j < 3; j++) {
403  triangle.getCorner(j, X1[0], X1[1]);
404 
405  Warp->computeDenom(X1, p_);
406  Warp->warpX(X1, X2, p_);
407  evolRMS += (x_pos[i * 3 + j] - X2[0]) * (x_pos[i * 3 + j] - X2[0]) +
408  (y_pos[i * 3 + j] - X2[1]) * (y_pos[i * 3 + j] - X2[1]);
409  x_pos[i * 3 + j] = X2[0];
410  y_pos[i * 3 + j] = X2[1];
411  }
412  }
413  evolRMS = evolRMS / nb_corners;
414 }
415 
vpImageFilter::getGradX
static void getGradX(const vpImage< unsigned char > &I, vpImage< double > &dIx)
Definition: vpImageFilter.cpp:545
vpTemplateTracker::G
vpColVector G
Definition: vpTemplateTracker.h:100
vpTemplateTrackerZNCCInverseCompositional::initHessienDesired
void initHessienDesired(const vpImage< unsigned char > &I)
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:77
vpColVector::resize
void resize(const unsigned int i, const bool flagNullify=true)
Definition: vpColVector.h:243
vpTemplateTracker::useInverse
bool useInverse
Definition: vpTemplateTracker.h:123
vpTemplateTracker::taillef
unsigned int taillef
Definition: vpTemplateTracker.h:108
vpImage::getValue
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1419
vpTemplateTrackerWarp
Definition: vpTemplateTrackerWarp.h:57
vpImageFilter::getGradY
static void getGradY(const vpImage< unsigned char > &I, vpImage< double > &dIy)
Definition: vpImageFilter.cpp:562
vpTemplateTrackerZNCCInverseCompositional::moydIrefdp
vpColVector moydIrefdp
Definition: vpTemplateTrackerZNCCInverseCompositional.h:66
vpTemplateTracker::ptTemplate
vpTemplateTrackerPoint * ptTemplate
Definition: vpTemplateTracker.h:69
vpColVector::t
vpRowVector t() const
Definition: vpColVector.cpp:695
vpTemplateTrackerZNCC
Definition: vpTemplateTrackerZNCC.h:64
vpTemplateTrackerZNCCInverseCompositional::trackNoPyr
void trackNoPyr(const vpImage< unsigned char > &I)
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:245
vpTemplateTracker::dp
vpColVector dp
Definition: vpTemplateTracker.h:128
vpTemplateTrackerPoint::dy
double dy
Definition: vpTemplateTrackerHeader.h:73
vpTemplateTrackerZNCCInverseCompositional::x_pos
std::vector< double > x_pos
Definition: vpTemplateTrackerZNCCInverseCompositional.h:63
vpImageFilter::filter
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, const bool convolve=false)
Definition: vpImageFilter.cpp:125
vpTemplateTrackerZNCCInverseCompositional::initCompInverse
void initCompInverse(const vpImage< unsigned char > &I)
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:51
vpTemplateTracker::diverge
bool diverge
Definition: vpTemplateTracker.h:120
vpTemplateTrackerZNCCInverseCompositional::computeEvalRMS
void computeEvalRMS(const vpColVector &p)
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:391
vpTemplateTrackerZNCCInverseCompositional::initPosEvalRMS
void initPosEvalRMS(const vpColVector &p)
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:369
vpTemplateTracker::X2
vpColVector X2
Definition: vpTemplateTracker.h:132
vpTemplateTracker::fgdG
double * fgdG
Definition: vpTemplateTracker.h:110
vpTemplateTrackerZone::getTriangle
void getTriangle(unsigned int i, vpTemplateTrackerTriangle &T) const
Definition: vpTemplateTrackerZone.cpp:353
vpTemplateTracker::fgG
double * fgG
Definition: vpTemplateTracker.h:109
vpTemplateTracker::blur
bool blur
Definition: vpTemplateTracker.h:105
vpTemplateTrackerTriangle::getCorner
vpColVector getCorner(unsigned int i) const
Definition: vpTemplateTrackerTriangle.h:115
vpTemplateTracker::HLMdesire
vpMatrix HLMdesire
Definition: vpTemplateTracker.h:96
vpTemplateTrackerWarp::dWarp
virtual void dWarp(const vpColVector &X1, const vpColVector &X2, const vpColVector &ParamM, vpMatrix &dW)=0
vpImageFilter::getGradYGauss2D
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
Definition: vpImageFilter.cpp:690
vpTemplateTrackerZone::getNbTriangle
unsigned int getNbTriangle() const
Definition: vpTemplateTrackerZone.h:99
vpTemplateTrackerWarp::warpX
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
vpTemplateTrackerWarp::getParamInverse
virtual void getParamInverse(const vpColVector &ParamM, vpColVector &ParamMinv) const =0
vpTemplateTracker::templateSize
unsigned int templateSize
Definition: vpTemplateTracker.h:72
vpTemplateTracker::X1
vpColVector X1
Definition: vpTemplateTracker.h:131
vpColVector
Implementation of column vector and the associated operations.
Definition: vpColVector.h:71
vpTemplateTracker::BI
vpImage< double > BI
Definition: vpTemplateTracker.h:136
vpTemplateTrackerZNCCInverseCompositional::y_pos
std::vector< double > y_pos
Definition: vpTemplateTrackerZNCCInverseCompositional.h:64
vpMatrix
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:103
vpTemplateTrackerWarp::getdW0
virtual void getdW0(const int &i, const int &j, const double &dy, const double &dx, double *dIdW)=0
vpTemplateTracker::nbIteration
unsigned int nbIteration
Definition: vpTemplateTracker.h:121
vpTemplateTracker::dW
vpMatrix dW
Definition: vpTemplateTracker.h:134
vpTemplateTracker::Hdesire
vpMatrix Hdesire
Definition: vpTemplateTracker.h:93
vpTemplateTrackerPoint::val
double val
Definition: vpTemplateTrackerHeader.h:74
vpTemplateTrackerWarp::pRondp
virtual void pRondp(const vpColVector &p1, const vpColVector &p2, vpColVector &pres) const =0
vpTemplateTracker::iterationMax
unsigned int iterationMax
Definition: vpTemplateTracker.h:116
vpTemplateTracker::HLMdesireInverse
vpMatrix HLMdesireInverse
Definition: vpTemplateTracker.h:98
vpTemplateTrackerZNCCInverseCompositional::threshold_RMS
double threshold_RMS
Definition: vpTemplateTrackerZNCCInverseCompositional.h:65
vpTemplateTracker::nbParam
unsigned int nbParam
Definition: vpTemplateTracker.h:114
vpImage::getHeight
unsigned int getHeight() const
Definition: vpImage.h:177
vpTemplateTracker::dIy
vpImage< double > dIy
Definition: vpTemplateTracker.h:138
vpMatrix::inverseByLU
vpMatrix inverseByLU() const
Definition: vpMatrix_lu.cpp:129
vpTemplateTrackerZNCCInverseCompositional::deletePosEvalRMS
void deletePosEvalRMS()
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:415
vpTemplateTrackerTriangle
Definition: vpTemplateTrackerTriangle.h:60
vpTemplateTrackerPoint::y
int y
Definition: vpTemplateTrackerHeader.h:72
vpTemplateTracker::dIx
vpImage< double > dIx
Definition: vpTemplateTracker.h:137
vpTemplateTracker::lambdaDep
double lambdaDep
Definition: vpTemplateTracker.h:115
vpTemplateTrackerPoint::dW
double * dW
Definition: vpTemplateTrackerHeader.h:75
vpTemplateTracker::p
vpColVector p
Definition: vpTemplateTracker.h:127
vpMatrix::computeHLM
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:5070
vpTemplateTracker::zoneTracked
vpTemplateTrackerZone * zoneTracked
Definition: vpTemplateTracker.h:87
vpTemplateTrackerZNCCInverseCompositional::vpTemplateTrackerZNCCInverseCompositional
vpTemplateTrackerZNCCInverseCompositional(vpTemplateTrackerWarp *warp)
Definition: vpTemplateTrackerZNCCInverseCompositional.cpp:44
vpImage< unsigned char >
vpTemplateTrackerZNCCInverseCompositional::evolRMS
double evolRMS
Definition: vpTemplateTrackerZNCCInverseCompositional.h:62
vpTemplateTrackerPoint::x
int x
Definition: vpTemplateTrackerHeader.h:72
vpTemplateTrackerPoint::dx
double dx
Definition: vpTemplateTrackerHeader.h:73
vpTemplateTracker::Warp
vpTemplateTrackerWarp * Warp
Definition: vpTemplateTracker.h:125
vpImage::getWidth
unsigned int getWidth() const
Definition: vpImage.h:238
vpImageFilter::getGradXGauss2D
static void getGradXGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIx, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
Definition: vpImageFilter.cpp:672
vpTemplateTrackerZNCCInverseCompositional::compoInitialised
bool compoInitialised
Definition: vpTemplateTrackerZNCCInverseCompositional.h:60