Visual Servoing Platform  version 3.2.0 under development (2019-01-22)
vpTemplateTrackerSSDESM.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/vpTemplateTrackerSSDESM.h>
42 
44  : vpTemplateTrackerSSD(warp), compoInitialised(false), HDir(), HInv(), HLMDir(), HLMInv(), GDir(), GInv()
45 {
46  useCompositionnal = false;
47  useInverse = false;
48 
49  if (!Warp->isESMcompatible())
50  std::cerr << "The selected warp function is not appropriate for the ESM "
51  "algorithm..."
52  << std::endl;
53 
60 }
61 
63 
65 {
66  // std::cout<<"Initialise precomputed value of ESM with templateSize: "<<
67  // templateSize<<std::endl;
69  int i, j;
70  // direct
71  for (unsigned int point = 0; point < templateSize; point++) {
72  i = ptTemplate[point].y;
73  j = ptTemplate[point].x;
74  ptTemplateCompo[point].dW = new double[2 * nbParam];
75  X1[0] = j;
76  X1[1] = i;
77  Warp->computeDenom(X1, p);
78  Warp->getdWdp0(i, j, ptTemplateCompo[point].dW);
79  }
80 
81  // inverse
82  HInv = 0;
83  for (unsigned int point = 0; point < templateSize; point++) {
84  i = ptTemplate[point].y;
85  j = ptTemplate[point].x;
86 
87  X1[0] = j;
88  X1[1] = i;
89  Warp->computeDenom(X1, p);
90  ptTemplate[point].dW = new double[nbParam];
91  Warp->getdW0(i, j, ptTemplate[point].dy, ptTemplate[point].dx, ptTemplate[point].dW);
92 
93  for (unsigned int it = 0; it < nbParam; it++)
94  for (unsigned int jt = 0; jt < nbParam; jt++)
95  HInv[it][jt] += ptTemplate[point].dW[it] * ptTemplate[point].dW[jt];
96  }
98 
99  compoInitialised = true;
100 }
101 
103 {
104  if (blur)
108 
109  double IW, dIWx, dIWy;
110  double Tij;
111  unsigned int iteration = 0;
112  int i, j;
113  double i2, j2;
114  double alpha = 2.;
115  do {
116  unsigned int Nbpoint = 0;
117  double erreur = 0;
118  dp = 0;
119  HDir = 0;
120  GDir = 0;
121  GInv = 0;
122  Warp->computeCoeff(p);
123  for (unsigned int point = 0; point < templateSize; point++) {
124  i = ptTemplate[point].y;
125  j = ptTemplate[point].x;
126  X1[0] = j;
127  X1[1] = i;
128 
129  Warp->computeDenom(X1, p);
130  Warp->warpX(X1, X2, p);
131 
132  j2 = X2[0];
133  i2 = X2[1];
134  if ((i2 >= 0) && (j2 >= 0) && (i2 < I.getHeight() - 1) && (j2 < I.getWidth() - 1)) {
135  // INVERSE
136  Tij = ptTemplate[point].val;
137  if (!blur)
138  IW = I.getValue(i2, j2);
139  else
140  IW = BI.getValue(i2, j2);
141  Nbpoint++;
142  double er = (Tij - IW);
143  for (unsigned int it = 0; it < nbParam; it++)
144  GInv[it] += er * ptTemplate[point].dW[it];
145 
146  erreur += er * er;
147 
148  // DIRECT
149  // dIWx=dIx.getValue(i2,j2);
150  // dIWy=dIy.getValue(i2,j2);
151 
152  dIWx = dIx.getValue(i2, j2) + ptTemplate[point].dx;
153  dIWy = dIy.getValue(i2, j2) + ptTemplate[point].dy;
154 
155  // Calcul du Hessien
156  // Warp->dWarp(X1,X2,p,dW);
157  Warp->dWarpCompo(X1, X2, p, ptTemplateCompo[point].dW, dW);
158 
159  double *tempt = new double[nbParam];
160  for (unsigned int it = 0; it < nbParam; it++)
161  tempt[it] = dW[0][it] * dIWx + dW[1][it] * dIWy;
162 
163  for (unsigned int it = 0; it < nbParam; it++)
164  for (unsigned int jt = 0; jt < nbParam; jt++)
165  HDir[it][jt] += tempt[it] * tempt[jt];
166 
167  for (unsigned int it = 0; it < nbParam; it++)
168  GDir[it] += er * tempt[it];
169  delete[] tempt;
170  }
171  }
172  if (Nbpoint == 0) {
173  // std::cout<<"plus de point dans template suivi"<<std::endl;
174  throw(vpTrackingException(vpTrackingException::notEnoughPointError, "No points in the template"));
175  }
176 
178 
179  try {
180  // dp=(HLMInv+HLMDir).inverseByLU()*(GInv+GDir);
181  // dp=HLMInv.inverseByLU()*GInv+HLMDir.inverseByLU()*GDir;
182  // dp=HLMInv.inverseByLU()*GInv;
183  dp = (HLMDir).inverseByLU() * (GDir);
184  } catch (const vpException &e) {
185  // std::cout<<"probleme inversion"<<std::endl;
186  throw(e);
187  }
188 
189  dp = gain * dp;
190  if (useBrent) {
191  alpha = 2.;
192  computeOptimalBrentGain(I, p, erreur / Nbpoint, dp, alpha);
193  dp = alpha * dp;
194  }
195 
196  // Warp->pRondp(p,dp,p);
197  p += dp;
198  iteration++;
199 
200  } while ((iteration < iterationMax));
201 
202  nbIteration = iteration;
203 }
204 
205 /*void vpTemplateTrackerSSDESM::InitCompInverse(vpImage<unsigned char> &I)
206 {
207  ptTempateCompo=new vpTemplateTrackerPointCompo[taille_template];
208  int i,j;
209  for(int point=0;point<taille_template;point++)
210  {
211  i=ptTemplate[point].y;
212  j=ptTemplate[point].x;
213  ptTempateCompo[point].dW=new double[2*nbParam];
214  Warp->getdWdp0(i,j,ptTempateCompo[point].dW);
215  }
216 
217 }
218 
219 void vpTemplateTrackerSSDESM::track(vpImage<unsigned char> &I)
220 {
221  double erreur=0,erreur_prec=1e38;
222  int Nbpoint=0;
223 
224  int taillefiltre=taille_filtre_dgaussien;
225  double *fg=new double[taillefiltre+1] ;
226  getGaussianDerivativeKernel(fg, taillefiltre) ;
227  getGradX(I, dIx,fg,taillefiltre);
228  getGradY(I, dIy,fg,taillefiltre);
229  delete[] fg;
230 
231  vpColVector dpinv(nbParam);
232  double lambda=lambdaDep;
233  double IW,dIWx,dIWy;
234  double Tij;
235  int iteration=0;
236  int i,j;
237  double i2,j2;
238  vpTemplateTrackerPoint *pt;
239  do
240  {
241  Nbpoint=0;
242  erreur_prec=erreur;
243  erreur=0;
244  dp=0;
245  HDir=0;
246  GDir=0;
247  GInv=0;
248  Warp->ComputeCoeff(p);
249  for(int point=0;point<taille_template;point++)
250  {
251  pt=&ptTemplate[point];
252  i=pt->y;
253  j=pt->x;
254  X1[0]=j;X1[1]=i;
255 
256  Warp->ComputeDenom(X1,p);
257  Warp->warpX(X1,X2,p);
258 
259  j2=X2[0];i2=X2[1];
260  if((j2<I.getWidth())&&(i2<I.getHeight())&&(i2>0)&&(j2>0))
261  {
262  //INVERSE
263  Tij=pt->val;
264  IW=I.getPixelBI(j2,i2);
265  Nbpoint++;
266  double er=(Tij-IW);
267  for(int it=0;it<nbParam;it++)
268  GInv[it]+=er*ptTemplate[point].dW[it];
269 
270  erreur+=er*er;
271 
272  //DIRECT COMPO
273  Tij=ptTemplate[point].val;
274  IW=I.getPixelBI(j2,i2);
275  dIWx=dIx.getPixelBI(j2,i2);
276  dIWy=dIy.getPixelBI(j2,i2);
277  Nbpoint++;
278  Warp->dWarpCompo(X1,X2,p,ptTempateCompo[point].dW,dW);
279  double *tempt=new double[nbParam];
280  for(int it=0;it<nbParam;it++)
281  tempt[it] =dW[0][it]*dIWx+dW[1][it]*dIWy;
282 
283  for(int it=0;it<nbParam;it++)
284  for(int jt=0;jt<nbParam;jt++)
285  HDir[it][jt]+=tempt[it]*tempt[jt];
286 
287  for(int it=0;it<nbParam;it++)
288  GDir[it]+=er*tempt[it];
289 
290  delete[] tempt;
291  }
292 
293 
294  }
295  if(Nbpoint==0)std::cout<<"plus de point dans template suivi"<<std::endl;
296  try
297  {
298  dp=(HInv+HDir).inverseByLU()*(GInv+GDir);
299  }
300  catch(...)
301  {
302  std::cout<<"probleme inversion"<<std::endl;
303  break;
304  }
305 
306  if(Compare)write_infos(p,erreur/Nbpoint);
307 
308  p+=Gain*dp;
309  iteration++;
310 
311  }
312  while( (iteration < IterationMax));
313 
314  NbIteration=iteration;
315 }*/
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
vpTemplateTrackerSSDESM(vpTemplateTrackerWarp *warp)
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
Definition: vpArray2D.h:171
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
void trackNoPyr(const vpImage< unsigned char > &I)
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.
void initHessienDesired(const vpImage< unsigned char > &I)
vpImage< double > dIx
vpImage< double > dIy
virtual bool isESMcompatible() const =0
virtual void getdW0(const int &i, const int &j, const double &dy, const double &dx, double *dIdW)=0
unsigned int nbIteration
vpTemplateTrackerPointCompo * ptTemplateCompo
void initCompInverse(const vpImage< unsigned char > &I)
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
void resize(const unsigned int i, const bool flagNullify=true)
Definition: vpColVector.h:244