Visual Servoing Platform  version 3.2.0 under development (2019-01-22)
vpTemplateTrackerSSDForwardCompositional.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 <visp3/core/vpImageFilter.h>
41 #include <visp3/tt/vpTemplateTrackerSSDForwardCompositional.h>
42 
44  : vpTemplateTrackerSSD(warp), compoInitialised(false)
45 {
46 }
47 
49 {
50  // std::cout<<"Initialise precomputed value of Compositionnal
51  // Direct"<<std::endl;
52  for (unsigned int point = 0; point < templateSize; point++) {
53  int i = ptTemplate[point].y;
54  int j = ptTemplate[point].x;
55  ptTemplate[point].dW = new double[2 * nbParam];
56  X1[0] = j;
57  X1[1] = i;
58  Warp->computeDenom(X1, p);
59  Warp->getdWdp0(i, j, ptTemplate[point].dW);
60  }
61  compoInitialised = true;
62 }
63 
65 
67 {
68  if (!compoInitialised)
69  std::cout << "Compositionnal tracking no initialised\nUse "
70  "InitCompo(vpImage<unsigned char> &I) function"
71  << std::endl;
72 
73  if (blur)
77 
78  dW = 0;
79 
80  double lambda = lambdaDep;
81  double IW, dIWx, dIWy;
82  double Tij;
83  unsigned int iteration = 0;
84  int i, j;
85  double i2, j2;
86  double alpha = 2.;
87  do {
88  unsigned int Nbpoint = 0;
89  double erreur = 0;
90  G = 0;
91  H = 0;
92  Warp->computeCoeff(p);
93  for (unsigned int point = 0; point < templateSize; point++) {
94  i = ptTemplate[point].y;
95  j = ptTemplate[point].x;
96  X1[0] = j;
97  X1[1] = i;
98 
99  Warp->computeDenom(X1, p);
100  Warp->warpX(X1, X2, p);
101 
102  j2 = X2[0];
103  i2 = X2[1];
104  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
105  Tij = ptTemplate[point].val;
106  if (!blur)
107  IW = I.getValue(i2, j2);
108  else
109  IW = BI.getValue(i2, j2);
110  dIWx = dIx.getValue(i2, j2);
111  dIWy = dIy.getValue(i2, j2);
112  Nbpoint++;
113  // Calcul du Hessien
114  /*Warp->dWarp(X1,X2,p,dW);
115  double *tempt=new double[nbParam];
116  for(int it=0;it<nbParam;it++)
117  tempt[it]=dW[0][it]*dIWx+dW[1][it]*dIWy;*/
118 
119  Warp->dWarpCompo(X1, X2, p, ptTemplate[point].dW, dW);
120 
121  double *tempt = new double[nbParam];
122  for (unsigned int it = 0; it < nbParam; it++)
123  tempt[it] = dW[0][it] * dIWx + dW[1][it] * dIWy;
124 
125  for (unsigned int it = 0; it < nbParam; it++)
126  for (unsigned int jt = 0; jt < nbParam; jt++)
127  H[it][jt] += tempt[it] * tempt[jt];
128 
129  double er = (Tij - IW);
130  for (unsigned int it = 0; it < nbParam; it++)
131  G[it] += er * tempt[it];
132 
133  erreur += (er * er);
134  delete[] tempt;
135  }
136  }
137  if (Nbpoint == 0) {
138  // std::cout<<"plus de point dans template suivi"<<std::endl;
139  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
140  }
141 
142  vpMatrix::computeHLM(H, lambda, HLM);
143 
144  try {
145  dp = 1. * HLM.inverseByLU() * G;
146  } catch (const vpException &e) {
147  // std::cout<<"probleme inversion"<<std::endl;
148  throw(e);
149  }
150 
151  dp = gain * dp;
152  if (useBrent) {
153  alpha = 2.;
154  computeOptimalBrentGain(I, p, erreur / Nbpoint, dp, alpha);
155  dp = alpha * dp;
156  }
157  Warp->pRondp(p, dp, p);
158  // p+=Gain*dp;
159  iteration++;
160  } while (/*( erreur_prec-erreur<50) &&*/ (iteration < iterationMax));
161 
162  nbIteration = iteration;
163 }
virtual void dWarpCompo(const vpColVector &X1, const vpColVector &X2, const vpColVector &ParamM, const double *dwdp0, vpMatrix &dW)=0
unsigned int getWidth() const
Definition: vpImage.h:239
vpTemplateTrackerPoint * ptTemplate
virtual void warpX(const int &i, const int &j, double &i2, double &j2, const vpColVector &ParamM)=0
void computeOptimalBrentGain(const vpImage< unsigned char > &I, vpColVector &tp, double tMI, vpColVector &direction, double &alpha)
error that can be emited by ViSP classes.
Definition: vpException.h:71
static void getGradYGauss2D(const vpImage< unsigned char > &I, vpImage< double > &dIy, const double *gaussianKernel, const double *gaussianDerivativeKernel, unsigned int size)
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
virtual void pRondp(const vpColVector &p1, const vpColVector &p2, vpColVector &pres) const =0
Type getValue(unsigned int i, unsigned int j) const
Definition: vpImage.h:1420
Error that can be emited by the vpTracker class and its derivates.
vpImage< double > dIx
vpImage< double > dIy
unsigned int nbIteration
vpMatrix inverseByLU() const
unsigned int getHeight() const
Definition: vpImage.h:178
static void filter(const vpImage< double > &I, vpImage< double > &Iu, vpImage< double > &Iv, const vpMatrix &M, const bool convolve=false)
virtual void getdWdp0(const int &i, const int &j, double *dIdW)=0
vpTemplateTrackerWarp * Warp
static void computeHLM(const vpMatrix &H, const double &alpha, vpMatrix &HLM)
Definition: vpMatrix.cpp:5071