Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TrackerBaseAffine.cpp
1 /*
2 This file is part of the BIAS library (Basic ImageAlgorithmS).
3 
4 Copyright (C) 2003-2009 (see file CONTACT for details)
5  Multimediale Systeme der Informationsverarbeitung
6  Institut fuer Informatik
7  Christian-Albrechts-Universitaet Kiel
8 
9 
10 BIAS is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 2.1 of the License, or
13 (at your option) any later version.
14 
15 BIAS is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Lesser General Public License for more details.
19 
20 You should have received a copy of the GNU Lesser General Public License
21 along with BIAS; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
24 
25 #include "TrackerBaseAffine.hh"
26 #include "MathAlgo/SVD.hh"
27 #include "Base/Image/ImageIO.hh"
28 #include "Base/Common/FileHandling.hh"
29 #include "Geometry/HMatrix.hh"
30 
31 using namespace BIAS;
32 using namespace std;
33 
34 const double DEFAULT_MAX_AFFINECHANGE = 0.3;
35 
36 
37 namespace BIAS{
38 
39 
40 void Compute4by1ErrorVector(KLT_TYPE* imgdiff,
41  KLT_TYPE* gradx,
42  KLT_TYPE* grady,
43  const int hw,
44  const int hh,
46 
47 void Compute4by4GradientMatrix(KLT_TYPE* gradx,
48  KLT_TYPE* grady,
49  const int hw,
50  const int hh,
52 
53 void Compute8by1ErrorVector(KLT_TYPE* imgdiff,
54  KLT_TYPE* gradx,
55  KLT_TYPE* grady,
56  KLT_TYPE a11,
57  KLT_TYPE a12,
58  KLT_TYPE a21,
59  KLT_TYPE a22,
60  KLT_TYPE tx,
61  KLT_TYPE ty,
62  KLT_TYPE h31,
63  KLT_TYPE h32,
64  const int hw,
65  const int hh,
67 
68 void Compute8by8GradientMatrix(KLT_TYPE* gradx,
69  KLT_TYPE* grady,
70  KLT_TYPE a11,
71  KLT_TYPE a12,
72  KLT_TYPE a21,
73  KLT_TYPE a22,
74  KLT_TYPE tx,
75  KLT_TYPE ty,
76  KLT_TYPE h31,
77  KLT_TYPE h32,
78  const int hw,
79  const int hh,
81 
84 
85 void Compute6by1ErrorVector(KLT_TYPE* imgdiff,
86  KLT_TYPE* gradx,
87  KLT_TYPE* grady,
88  const int hw,
89  const int hh,
91 
92 void Compute6by6GradientMatrix(KLT_TYPE* gradx,
93  KLT_TYPE* grady,
94  const int hw,
95  const int hh,
97 
98 
99 template <class StorageType>
101 TrackerBaseAffine(bool only4params, bool useHomography)
102 {
103  _MaxAffineChange = DEFAULT_MAX_AFFINECHANGE;
104  _SimilarityTransformOnly = only4params;
105  UseHomography_ = useHomography;
106 }
107 
108 template <class StorageType>
111  Vector2<KLT_TYPE>& result, KLT_TYPE& error,
112  int &iter,
113  const Matrix2x2<KLT_TYPE>& AffinePred,
114  Matrix2x2<KLT_TYPE>& AffineResult)
115 {
116 
117 #ifdef BIAS_DEBUG
118  if (AffinePred.NormFrobenius()<1e-10) {
119  BIASERR("affine guess seems unitialized, set at least identity");
120  }
121 #endif
122  residuumMAD_ = 0;
123  residuumMSD_ = 0;
124 
125  int res = TrackAffine_(p1.GetData(), p2.GetData(), AffinePred,
126  result.GetData(), AffineResult,
127  error, iter, Cov_, residuumMAD_, residuumMSD_,
128  _SimilarityTransformOnly, true);
129  return res;
130 }
131 
132 // ***************** parts, which are not under LGPL *************************
133 // The function TrackAffine_ is based on the affine tracker by Thorsten
134 // Thormaehlen. Some of the worker functions have been rewritten, some
135 // have been adapted, see below for the conditions on this software.
136 
137 template <class StorageType>
139 TrackAffine_(KLT_TYPE p1[2], KLT_TYPE p2[2],
140  const Matrix2x2<KLT_TYPE>& p2A,
141  KLT_TYPE result[2], Matrix2x2<KLT_TYPE>& resultA,
142  KLT_TYPE& error, int &iteration, Matrix<KLT_TYPE>& Cov,
143  KLT_TYPE& mad, KLT_TYPE& msd,
144  bool SimilarityTransformOnly, bool ComputeCovariance)
145 
146 {
147 
148  if (_HalfWinSize != _LastHalfWinSize){
149  this->Resize_(_HalfWinSize);
150  }
151 #ifdef BIAS_DEBUG
152  static int featurecounter = 0;
153  featurecounter++;
154 #endif
155 
156 
157  //this->AddDebugLevel(D_TRACKERB_AFFINE | D_TRACKERB_KLT);
158  BIASDOUT(D_TRACKERB_AFFINE, "Called Affine for "
159  <<p1[0]<<";"<<p1[1]<<", guess is "<<p2[0]<<";"<<p2[1]
160  <<" affine guess is "<< p2A);
161 
162  const int width = _WinSize;
163  const int height = _WinSize;
164  const int hw = width/2;
165  const int hh = height/2;
166 
167  const KLT_TYPE fhw = hw;
168  const KLT_TYPE fhh = hh;
169 
170  if (_minx1 >= p1[0] || _miny1 >= p1[1] || _maxx1 < p1[0] || _maxy1 < p1[1]){
171  BIASDOUT(D_TRACKERB_KLT, "original point is too close to border: "
172  <<p1[0]<<";"<<p1[1]);
173  return TRACKER_OUTOFIMAGE;
174  }
175 
176 
177  const int max_iterations = _MaxIterations;
178  const KLT_TYPE maxerr = _MaxError;
179  const KLT_TYPE th_aff = _MaxAffineChange;
180 
181  KLT_TYPE dx=0, dy=0;
182  iteration = 0;
183  int status = 0;
184 
185  Matrix<KLT_TYPE> T(UseHomography_?8:SimilarityTransformOnly?4:6,
186  UseHomography_?8:SimilarityTransformOnly?4:6);
187  Vector<KLT_TYPE> a(UseHomography_?8:SimilarityTransformOnly?4:6);
188 
189  Cov.newsize(SimilarityTransformOnly?4:6,SimilarityTransformOnly?4:6);
190 
191 
192  resultA = p2A;
193  KLT_TYPE *Axx = &(resultA[0][0]);
194  KLT_TYPE *Axy = &(resultA[0][1]);
195  KLT_TYPE *Ayx = &(resultA[1][0]);
196  KLT_TYPE *Ayy = &(resultA[1][1]);
197 
198  result[0] = p2[0];
199  result[1] = p2[1];
200 
201  KLT_TYPE h31=0,h32=0;
202 
203  bool convergence = false;
204 
205 
206  // Allocate memory for windows
207  KLT_TYPE *imgdiff = new KLT_TYPE[width * height];
208  KLT_TYPE *gradx = new KLT_TYPE[width * height];
209  KLT_TYPE *grady = new KLT_TYPE[width * height];
210 
211  // compute reference window only ONCE
212  BilinearRegion1_(p1[0], p1[1], hh);
213 
214  do {
215  KLT_TYPE ul_x = (result[0] + *Axy * fhh - *Axx * fhw)/(-h31*fhw+h32*fhh+1);
216  KLT_TYPE ul_y = (result[1] + *Ayy * fhh - *Ayx * fhw)/(-h31*fhw+h32*fhh+1);
217  KLT_TYPE ll_x = (result[0] - *Axx * fhw - *Axy * fhh)/(-h31*fhw-h32*fhh+1);
218  KLT_TYPE ll_y = (result[1] - *Ayx * fhw - *Ayy * fhh)/(-h31*fhw-h32*fhh+1);
219  KLT_TYPE ur_x = (result[0] + *Axx * fhw + *Axy * fhh)/(h31*fhw+h32*fhh+1);
220  KLT_TYPE ur_y = (result[1] + *Ayx * fhw + *Ayy * fhh)/(h31*fhw+h32*fhh+1);
221  KLT_TYPE lr_x = (result[0] + *Axx * fhw - *Axy * fhh)/(h31*fhw-h32*fhh+1);
222  KLT_TYPE lr_y = (result[1] + *Ayx * fhw - *Ayy * fhh)/(h31*fhw-h32*fhh+1);
223 
224  if ( ul_x < _minx2 || ul_x > _maxx2 ||
225  ll_x < _minx2 || ll_x > _maxx2 ||
226  ur_x < _minx2 || ur_x > _maxx2 ||
227  lr_x < _minx2 || lr_x > _maxx2 ||
228  ul_y < _miny2 || ul_y > _maxy2 ||
229  ll_y < _miny2 || ll_y > _maxy2 ||
230  ur_y < _miny2 || ur_y > _maxy2 ||
231  lr_y < _miny2 || lr_y > _maxy2) {
232  status = TRACKER_OUTOFIMAGE;
233  BIASDOUT( D_TRACKERB_KLT, "out of image");
234  break;
235  }
236 
237  // compute deltaI and also brightness normalization factor dev2_
238  ComputeIntensityDifferenceAffine_(result[0], result[1],
239  *Axx, *Ayx , *Axy, *Ayy, h31, h32,
240  hw, hh, imgdiff);
241 
242  // compute gradients and apply dev2_ to gradients
243  GetGradientWinAffine_(result[0], result[1],
244  *Axx, *Ayx , *Axy, *Ayy, h31, h32,
245  hw, hh, gradx, grady);
246 
247 
248 
249 
250 
251 #ifdef BIAS_DEBUG
252  if (this->GetDebugLevel() & D_TRACKERB_AFFINE) {
253  stringstream ss;
254  ss<<"Aff_"<<setw(3)<<setfill('0')<<int(featurecounter)<<"-"
255  <<setw(4)<<setfill('0')<<int(iteration)<<"-";
256  string imname = ss.str();
257  Image<float> im1(2*hw+1, 2*hh+1), im2(2*hw+1, 2*hh+1);
258  im1.SetZero();
259  im2.SetZero();
260  for (int y=0; y<2*hh+1; y++) {
261  for (int x=0; x<2*hw+1; x++) {
262  im1.GetImageDataArray()[y][x] = float(_bl1[y][x]);
263  im2.GetImageDataArray()[y][x] = float(_bl2[y][x]);
264  }
265  }
266  //ImageIO::Save(imname+FileHandling::toString("1"), im1);
267  //ImageIO::Save(imname+FileHandling::toString("2"), im2);
268  ImageIO::Save(imname+FileHandling::toString("1"), im1);
269  ImageIO::Save(imname+FileHandling::toString("2"), im2);
270  }
271 #endif
272 
273 
274  if (UseHomography_) {
275  Compute8by1ErrorVector(imgdiff, gradx, grady,
276  *Axx, *Axy, *Ayx, *Ayy,
277  //result[0]-p1[0], result[1]-p1[1], // correct ?
278  0,0,
279  h31, h32,
280  hw, hh, a);
281  Compute8by8GradientMatrix(gradx, grady,
282  *Axx, *Axy, *Ayx, *Ayy,
283  0,0,
284  //result[0]-p1[0], result[1]-p1[1], // correct ?
285  h31, h32,
286  hw, hh, T);
287 
288  SVD testsvd(T);
289  //cout<<"singular values for T: "<<testsvd.GetS()<<endl;
290  Vector<double> par(8);
291  status = testsvd.Solve(T,a,par);
292  a = par;
293  //status = GaussJordanElimination(T,a);
294 
295 
296 
297 
298 
299 
300  if (status==TRACKER_NOSTRUCTURE) {
301  BIASERR("structure degeneracy for homography estimation");
302  }
303  //cout<<"update for params is "<<a<<endl;
304  *Axx += a[0];
305  *Ayx += a[1];
306  *Axy += a[2];
307  *Ayy += a[3];
308 
309  result[0] += dx = a[4];
310  result[1] += dy = a[5];
311 
312  h31 += a[6];
313  h32 += a[7];
314  } else {
315 
316 
317 
318 
319  // // weight center stronger than outer regions !
320 // const double sigmainv = 1.0;
321 // double* pimgdiff = imgdiff;
322 // double* pgradx = gradx;
323 // double* pgrady = grady;
324 // for (int wy=-hh; wy<hh; wy++) {
325 // for (int wx=-hw; wx<hw; wx++) {
326 // const double weightfactor = exp(-sigmainv*double(wx*wx+wy*wy));
327 // *pimgdiff *= weightfactor;
328 // pimgdiff++;
329 
330 // *pgradx *= weightfactor;
331 // pgradx++;
332 
333 // *pgrady *= weightfactor;
334 // pgrady++;
335 // }
336 // }
337 
338 
339 
340 
341  if (SimilarityTransformOnly){
342  Compute4by1ErrorVector(imgdiff, gradx, grady, hw, hh, a);
343  Compute4by4GradientMatrix(gradx, grady, hw, hh, T);
344 
345  status = GaussJordanElimination(T,a);
346 
347  *Axx += a[0];
348  *Ayx += a[1];
349  *Ayy = *Axx;
350  *Axy = -(*Ayx);
351 
352  result[0] += dx = a[2];
353  result[1] += dy = a[3];
354  } else {
355  Compute6by1ErrorVector(imgdiff, gradx, grady, hw, hh, a);
356  Compute6by6GradientMatrix(gradx, grady, hw, hh, T);
357 
358  status = GaussJordanElimination(T,a);
359 
360  *Axx += a[0];
361  *Ayx += a[1];
362  *Axy += a[2];
363  *Ayy += a[3];
364 
365  result[0] += dx = a[4];
366  result[1] += dy = a[5];
367  }
368  }
369 
370  /*
371  double posdiag =
372  (*Axx + *Axy) * (*Axx + *Axy) + (*Ayx + *Ayy) * (*Ayx + *Ayy);
373 
374  double negdiag =
375  (*Axy - *Axx) * (*Axy - *Axx) + (*Ayy - *Ayx) * (*Ayy - *Ayx);
376 
377  double ratio = (posdiag>negdiag) ? negdiag / posdiag : posdiag / negdiag;
378 
379  // factor 3 is allowed
380  if (ratio < 0.09) {
381  status = TRACKER_ODDAFFINEWARP;
382  BIASDOUT(D_TRACKERB_KLT, "diag ratio is too big: "<< sqrt(ratio));
383  break;
384  }
385 
386  double scalex = (*Axx * *Axx) + (*Ayx * *Ayx);
387  double scaley = (*Ayy * *Ayy) + (*Axy * *Axy);
388 
389 
390  ratio = (scalex>scaley) ? scaley / scalex : scalex / scaley;
391  // factor 3 is allowed
392  if (ratio < 0.09) {
393  status = TRACKER_ODDAFFINEWARP;
394  BIASDOUT(D_TRACKERB_KLT, "axes ratio is too big: "<< sqrt(ratio));
395  break;
396  }
397  */
398 
399  // old corner - new corner, how far have we moved this iteration
400  ul_x -= (result[0] + *Axy * fhh - *Axx * fhw)/(-h31*fhw+h32*fhh+1);
401  ul_y -= (result[1] + *Ayy * fhh - *Ayx * fhw)/(-h31*fhw+h32*fhh+1);
402  ll_x -= (result[0] - *Axx * fhw - *Axy * fhh)/(-h31*fhw-h32*fhh+1);
403  ll_y -= (result[1] - *Ayx * fhw - *Ayy * fhh)/(-h31*fhw-h32*fhh+1);
404  ur_x -= (result[0] + *Axx * fhw + *Axy * fhh)/(h31*fhw+h32*fhh+1);
405  ur_y -= (result[1] + *Ayx * fhw + *Ayy * fhh)/(h31*fhw+h32*fhh+1);
406  lr_x -= (result[0] + *Axx * fhw - *Axy * fhh)/(h31*fhw-h32*fhh+1);
407  lr_y -= (result[1] + *Ayx * fhw - *Ayy * fhh)/(h31*fhw-h32*fhh+1);
408 
409 
410  convergence = (fabs(dx) < maxerr && fabs(dy) < maxerr &&
411  fabs(ul_x) < th_aff && fabs(ul_y) < th_aff &&
412  fabs(ll_x) < th_aff && fabs(ll_y) < th_aff &&
413  fabs(ur_x) < th_aff && fabs(ur_y) < th_aff &&
414  fabs(lr_x) < th_aff && fabs(lr_y) < th_aff);
415 
416  if (status == TRACKER_NOSTRUCTURE) {
417  BIASDOUT( D_TRACKERB_KLT, "no structure for tracking");
418  break;
419  }
420  iteration++;
421 
422  } while ( !convergence && iteration < max_iterations);
423  if (iteration >= max_iterations) {
424  status = TRACKER_TOOMANYITER;
425  BIASDOUT( D_TRACKERB_KLT, "too many iter: "<<iteration
426  <<" maxiter="<<max_iterations);
427  }
428  mad = 0;
429  msd = 0;
430  if ((status==TRACKER_TOOMANYITER) || (status==TRACKER_SUCCESS)) {
431 
432  /** if we compute the difference here, we must check the range again
433  ComputeIntensityDifferenceAffine(*_im1, *_im2, p1[0], p1[1],
434  result[0], result[1],
435  *Axx, *Ayx , *Axy, *Ayy,
436  hw, hh, imgdiff);
437  */
438  if (UseHomography_) {
439  HMatrix H;
440  H[0][0] = *Axx;
441  H[0][1] = *Axy;
442  H[1][0] = *Ayx;
443  H[1][1] = *Ayy;
444  H[0][2] = result[0];
445  H[1][2] = result[1];
446  H[2][0] = h31;
447  H[2][1] = h32;
448  H[2][2] = 1.0;
449  Matrix2x2<double> Jac;
450  H.GetJacobian(HomgPoint2D(0, 0), Jac); // correct ?
451  // cout<<Jac<<endl;
452  *Axx = Jac[0][0];
453  *Axy = Jac[0][1];
454  *Ayx = Jac[1][0];
455  *Ayy = Jac[1][1];
456  }
457 
458  const unsigned int wh = width*height;
459  register KLT_TYPE* fw = imgdiff;
460  for (unsigned int h=0; h<wh ; h++) {
461  const KLT_TYPE cur = (KLT_TYPE)fabsf((float) (*fw) );
462  fw++;
463  mad += cur;
464  msd += cur*cur;
465  }
466 
467  mad /= double(wh)*dev1_;
468  msd /= dev1_*dev1_;
469  double sigma = msd / (double)(wh - SimilarityTransformOnly?4:6);
470  msd /= double(wh);
471 
472  if (ComputeCovariance) {
473  if (UseHomography_)
474  T.newsize(SimilarityTransformOnly?4:6,SimilarityTransformOnly?4:6);
475  if (SimilarityTransformOnly)
476  Compute4by4GradientMatrix(gradx, grady, hw, hh, T);
477  else
478  Compute6by6GradientMatrix(gradx, grady, hw, hh, T);
479  SVD covsvd(T);
480  Cov = covsvd.Invert() * sigma;
481  } else {
482  // pseudo covariance (last step size)
483  Cov.SetZero();
484  for (int i=0; i<Cov.num_cols()-2; i++) {
485  Cov[i][i] = 0.01 * (dx * dx + dy * dy);
486  }
487  if (SimilarityTransformOnly) {
488  Cov[2][2] = dx * dx;
489  Cov[3][3] = dy * dy;
490  } else {
491  Cov[4][4] = dx * dx;
492  Cov[5][5] = dy * dy;
493  }
494  }
495  }
496 
497  error = (fabs(dx)>fabs(dy))?fabs(dx):fabs(dy);
498 
499  // Check whether window is out of bounds
500  if (result[0]<_minx2 || result[0] > _maxx2 ||
501  result[1]<_miny2 || result[1] > _maxy2)
502  status = TRACKER_OUTOFIMAGE;
503 
504 
505  // Free memory
506  delete imgdiff;
507  delete gradx;
508  delete grady;
509 
510  BIASDOUT( D_TRACKERB_KLT, "status at end is "<<status);
511 
512  // Return appropriate value
513  return status;
514 }
515 
516 template <class StorageType>
518 EvaluateResult_(KLT_TYPE& mad, KLT_TYPE& msd, Matrix<KLT_TYPE>& CovMatrix) {
519  // these values have already been computed in the track_ function !!!
520  if (_SimilarityTransformOnly) {
521  CovMatrix.newsize(4,4);
522  } else {
523  CovMatrix.newsize(6,6);
524  }
525  CovMatrix = Cov_;
526 
527  mad = residuumMAD_;
528  msd = residuumMSD_;
529 }
530 
531 
532 /**********************************************************************
533 * CONSISTENCY CHECK OF FEATURES BY AFFINE MAPPING (BEGIN)
534 *
535 * Created by: Thorsten Thormaehlen (University of Hannover) June 2004
536 * thormae@tnt.uni-hannover.de
537 *
538 * Permission is granted to any individual or institution to use, copy, modify,
539 * and distribute this part of the software, provided that this complete
540 * authorship and permission notice is maintained, intact, in all copies.
541 *
542 * This software is provided "as is" without express or implied warranty.
543 *
544 *
545 */
546 
548  Vector<KLT_TYPE>& b) {
549  const int n = a.GetRows();
550  Vector<int> indxc(n),indxr(n),ipiv(n);
551  int i,j,k,l,ll;
552  KLT_TYPE big,dum,pivinv,tmp;
553  int col = 0;
554  int row = 0;
555 
556  for (j=0;j<n;j++) ipiv[j]=0;
557  for (i=0;i<n;i++) {
558  big=0.0;
559  for (j=0;j<n;j++)
560  if (ipiv[j] != 1)
561  for (k=0;k<n;k++) {
562  if (ipiv[k] == 0) {
563  if (fabs(a[j][k]) >= big) {
564  big=fabs(a[j][k]);
565  row=j;
566  col=k;
567  }
568  } else if (ipiv[k] > 1) return TRACKER_NOSTRUCTURE;
569  }
570  ++(ipiv[col]);
571  if (row != col) {
572  for (l=0;l<n;l++) {
573  tmp = a[row][l];
574  a[row][l] = a[col][l];
575  a[col][l] = tmp;
576  }
577  tmp = b[row];
578  b[row] = b[col];
579  b[col] = tmp;
580  }
581  indxr[i]=row;
582  indxc[i]=col;
583  if (a[col][col] == 0.0) return TRACKER_NOSTRUCTURE;
584  pivinv=1.0/a[col][col];
585  a[col][col]=1.0;
586  for (l=0;l<n;l++) a[col][l] *= pivinv;
587  b[col] *= pivinv;
588  for (ll=0;ll<n;ll++)
589  if (ll != col) {
590  dum=a[ll][col];
591  a[ll][col]=0.0;
592  for (l=0;l<n;l++) a[ll][l] -= a[col][l]*dum;
593  b[ll] -= b[col]*dum;
594  }
595  }
596  for (l=n-1;l>=0;l--) {
597  if (indxr[l] != indxc[l])
598  for (k=0;k<n;k++) {
599  tmp = a[k][indxr[l]];
600  a[k][indxr[l]] = a[k][indxc[l]];
601  a[k][indxc[l]] = tmp;
602  }
603  }
604 
605  return TRACKER_SUCCESS;
606 }
607 
608 template <class StorageType>
610 GetGradientWinAffine_(const KLT_TYPE& x, const KLT_TYPE& y,
611  const KLT_TYPE& Axx, const KLT_TYPE& Ayx,
612  const KLT_TYPE& Axy, const KLT_TYPE& Ayy,
613  const KLT_TYPE& h31, const KLT_TYPE& h32,
614  const int hw, const int hh,
615  KLT_TYPE* out_gradx,
616  KLT_TYPE* out_grady) {
617  if (_ComputeFilteredImages) {
618 
619  double dev2 = dev2_;
620 
621  BIASWARNONCE("Inefficient code: affine sparse tracking, needs rewrite."
622  <<" either use precomputed gradients or implment this ");
623 
624 
625  // for some reason this approximation does not work (traking interrupts)
626  // dont debug this because it is very inefficient either !!!
627  BIASABORT;
628 
629 
630 
631 
632  bool AffineBrightnessInvariance = _AffineBrightnessInvariance;
633  _AffineBrightnessInvariance = false;
634 
635  register int i, j;
636  KLT_TYPE mi, mj;
637  Matrix<KLT_TYPE> I(2*hh+1, 2*hw+1, MatrixZero),
638  gx(2*hh+1, 2*hw+1, MatrixZero), gy(2*hh+1, 2*hw+1, MatrixZero);
639  for (j = -hh ; j <= hh ; j++) {
640  for (i = -hw ; i <= hw ; i++) {
641 
642 
643 
644  // FIXME: reeeeally inefficient!
645 
646 
647 
648  mi = (x + Axx * i + Axy * j) / (1.0 + h31*i + h32*j);
649  mj = (y + Ayx * i + Ayy * j) / (1.0 + h31*i + h32*j);
650 
651  BilinearRegion2_(mi, mj, 0);
652  gx[j+hh][i+hw] = _gx2[0][0];
653  gy[j+hh][i+hw] = _gy2[0][0];
654  I[j+hh][i+hw] = _bl2[0][0];
655  }
656  }
657  _AffineBrightnessInvariance = AffineBrightnessInvariance;
658  if (_AffineBrightnessInvariance) {
659  this->NormalizeRegion_(I, gx, gy);
660  }
661  dev2_ = dev2;
662  KLT_TYPE *pgx = gx.GetData(), *pgy = gy.GetData();
663  for (int c = (2*hh+1)*(2*hw+1); c>0; c--) {
664  *out_gradx++ = *pgx++;
665  *out_grady++ = *pgy++;
666  }
667  } else {
668  register int i, j;
669  KLT_TYPE mi, mj;
670  if (_AffineBrightnessInvariance) {
671  for (j = -hh ; j <= hh ; j++) {
672  for (i = -hw ; i <= hw ; i++) {
673  mi = (x + Axx * i + Axy * j) / (1.0 + h31*i + h32*j);
674  mj = (y + Ayx * i + Ayy * j) / (1.0 + h31*i + h32*j);
675 
676  // normalize by precomputed dev2_
677  *out_gradx++ = _gradim2x->FastBilinearInterpolationGrey(mi,
678  mj)*dev2_;
679  *out_grady++ = _gradim2y->FastBilinearInterpolationGrey(mi,
680  mj)*dev2_;
681  }
682  }
683  } else {
684  for (j = -hh ; j <= hh ; j++) {
685  for (i = -hw ; i <= hw ; i++) {
686  mi = (x + Axx * i + Axy * j) / (1.0 + h31*i + h32*j);
687  mj = (y + Ayx * i + Ayy * j) / (1.0 + h31*i + h32*j);
688  *out_gradx++ = _gradim2x->FastBilinearInterpolationGrey(mi, mj);
689  *out_grady++ = _gradim2y->FastBilinearInterpolationGrey(mi, mj);
690  }
691  }
692  }
693  }
694 }
695 
696 template <class StorageType>
698 ComputeIntensityDifferenceAffine_(const KLT_TYPE& x2, const KLT_TYPE& y2,
699  const KLT_TYPE& Axx, const KLT_TYPE& Ayx,
700  const KLT_TYPE& Axy, const KLT_TYPE& Ayy,
701  const KLT_TYPE& h31, const KLT_TYPE& h32,
702  const int hw, const int hh,
703  KLT_TYPE* imgdiff) {
704  register int i, j;
705  BIASASSERT(hh==hw);
706  KLT_TYPE* imgdifforig = imgdiff;
707  if (_ComputeFilteredImages) {
708  BIASERR("FIXME: reeeeally inefficient to interpolate always with size 1 !!!");
709  BIASABORT;
710  for (j = -hh ; j <= hh ; j++) {
711  for (i = -hw ; i <= hw ; i++) {
712  // FIXME: reeeeally inefficient to interpolate always with size 1 !!!
713  /// write specialized affine function
714  KLT_TYPE mi = (x2 + Axx * i + Axy * j) / (1.0 + h31*i + h32*j);
715  KLT_TYPE mj = (y2 + Ayx * i + Ayy * j) / (1.0 + h31*i + h32*j);
716  BilinearRegion2_(mi, mj, 0);
717  *imgdiff++ = _bl1[j+hh][i+hw] - _bl2[0][0];
718  }
719  }
720  } else {
721  mean1_ = mean2_ = 0;
722  BIASASSERT(_bl1.num_rows()==_bl1.num_cols());
723  BIASASSERT(_bl1.num_rows()==2*hw+1);
724  KLT_TYPE mi,mj;
725  KLT_TYPE *pbl1 = _bl1.GetData();
726  KLT_TYPE *pbl2 = _bl2.GetData();
727  for (j = -hh ; j <= hh ; j++) {
728  for (i = -hw ; i <= hw ; i++) {
729  mi = (x2 + Axx * i + Axy * j) / (1.0 + h31*i + h32*j);
730  mj = (y2 + Ayx * i + Ayy * j) / (1.0 + h31*i + h32*j);
731  mean2_ += *pbl2 = _im2->FastBilinearInterpolationGrey(mi, mj);
732  *imgdiff++ = *pbl1++ - *pbl2++;
733  }
734  }
735  if (_AffineBrightnessInvariance) {
736  pbl2 = _bl2.GetData();
737  const KLT_TYPE thearea = (2*hh+1)*(2*hw+1);
738  mean2_ /= thearea;
739  dev2_ = 0;
740  double diff;
741  for (j = -hh ; j <= hh ; j++)
742  for (i = -hw ; i <= hw ; i++) {
743  diff = (*pbl2++ -= mean2_);
744  dev2_ += diff * diff;
745  }
746  if (fabs(dev2_)>1e-10) dev2_ = sqrt(thearea / dev2_);
747  else dev2_ = 1.0;
748  imgdiff = imgdifforig;
749  pbl1 = _bl1.GetData();
750  pbl2 = _bl2.GetData();
751  for (j = -hh ; j <= hh ; j++) {
752  for (i = -hw ; i <= hw ; i++) {
753  *imgdiff++ =(*pbl1++) - (*pbl2++ *= dev2_);
754  }
755  }
756  }
757  }
758 }
759 
760 
761 void Compute6by6GradientMatrix(KLT_TYPE* gradx,
762  KLT_TYPE* grady,
763  const int hw,
764  const int hh,
765  Matrix<KLT_TYPE>& T) {
766  register int i, j;
767  KLT_TYPE gx, gy, gxx, gxy, gyy, x, y, xx, xy, yy;
768 
769 
770  // init with zero
771  for (j = 0 ; j < 6 ; j++) {
772  for (i = j ; i < 6 ; i++) {
773  T[j][i] = 0.0;
774  }
775  }
776 
777  for (j = -hh ; j <= hh ; j++) {
778  for (i = -hw ; i <= hw ; i++) {
779  gx = *gradx++;
780  gy = *grady++;
781  gxx = gx * gx;
782  gxy = gx * gy;
783  gyy = gy * gy;
784  x = i;
785  y = j;
786  xx = x * x;
787  xy = x * y;
788  yy = y * y;
789 
790  T[0][0] += xx * gxx;
791  T[0][1] += xx * gxy;
792  T[0][2] += xy * gxx;
793  T[0][3] += xy * gxy;
794  T[0][4] += x * gxx;
795  T[0][5] += x * gxy;
796 
797  T[1][1] += xx * gyy;
798  T[1][2] += xy * gxy;
799  T[1][3] += xy * gyy;
800  T[1][4] += x * gxy;
801  T[1][5] += x * gyy;
802 
803  T[2][2] += yy * gxx;
804  T[2][3] += yy * gxy;
805  T[2][4] += y * gxx;
806  T[2][5] += y * gxy;
807 
808  T[3][3] += yy * gyy;
809  T[3][4] += y * gxy;
810  T[3][5] += y * gyy;
811 
812  T[4][4] += gxx;
813  T[4][5] += gxy;
814 
815  T[5][5] += gyy;
816 
817 
818  }
819  }
820 
821  for (j = 0 ; j < 5 ; j++) {
822  for (i = j+1 ; i < 6 ; i++) {
823  T[i][j] = T[j][i];
824  }
825  }
826  // regularisation:
827  // double mean = 0;
828 // for (j = 0 ; j < 6 ; j++) {
829 // mean += fabs(T[j][j]);
830 // }
831 // // make sure we get no worse conditioning than 10^4
832 // mean *= 0.0001;
833 // for (j = 0 ; j < 6 ; j++) {
834 // T[j][j] += mean;
835 // }
836 
837 }
838 
839 void Compute6by1ErrorVector(KLT_TYPE* imgdiff,
840  KLT_TYPE* gradx,
841  KLT_TYPE* grady,
842  const int hw,
843  const int hh,
844  Vector<KLT_TYPE>& e) {
845  register int i, j;
846  register KLT_TYPE diff, diffgradx, diffgrady;
847 
848  /* Set values to zero */
849  for(i = 0; i < 6; i++) e[i] = 0.0;
850 
851  /* Compute values */
852  for (j = -hh ; j <= hh ; j++) {
853  for (i = -hw ; i <= hw ; i++) {
854  diff = *imgdiff++;
855  diffgradx = diff * (*gradx++);
856  diffgrady = diff * (*grady++);
857  e[0] += diffgradx * i;
858  e[1] += diffgrady * i;
859  e[2] += diffgradx * j;
860  e[3] += diffgrady * j;
861  e[4] += diffgradx;
862  e[5] += diffgrady;
863 
864  }
865  }
866 
867  for(i = 0; i < 6; i++) e[i] *= 0.5;
868 
869 }
870 
871 
872 void Compute4by4GradientMatrix(KLT_TYPE* gradx,
873  KLT_TYPE* grady,
874  const int hw,
875  const int hh,
876  Matrix<KLT_TYPE>& T)
877 {
878  register int i, j;
879  KLT_TYPE gx, gy, x, y;
880 
881 
882  // Set values to zero
883  for (j = 0 ; j < 4 ; j++) {
884  for (i = 0 ; i < 4 ; i++) {
885  T[j][i] = 0.0;
886  }
887  }
888 
889  for (j = -hh ; j <= hh ; j++) {
890  for (i = -hw ; i <= hw ; i++) {
891  gx = *gradx++;
892  gy = *grady++;
893  x = i;
894  y = j;
895  T[0][0] += (x*gx+y*gy) * (x*gx+y*gy);
896  T[0][1] += (x*gx+y*gy)*(x*gy-y*gx);
897  T[0][2] += (x*gx+y*gy)*gx;
898  T[0][3] += (x*gx+y*gy)*gy;
899 
900  T[1][1] += (x*gy-y*gx) * (x*gy-y*gx);
901  T[1][2] += (x*gy-y*gx)*gx;
902  T[1][3] += (x*gy-y*gx)*gy;
903 
904  T[2][2] += gx*gx;
905  T[2][3] += gx*gy;
906 
907  T[3][3] += gy*gy;
908  }
909  }
910 
911  for (j = 0 ; j < 3 ; j++) {
912  for (i = j+1 ; i < 4 ; i++) {
913  T[i][j] = T[j][i];
914  }
915  }
916 
917 }
918 
919 void Compute4by1ErrorVector(KLT_TYPE* imgdiff,
920  KLT_TYPE* gradx,
921  KLT_TYPE* grady,
922  const int hw,
923  const int hh,
924  Vector<KLT_TYPE>& e){
925  register int i, j;
926  register KLT_TYPE diff, diffgradx, diffgrady;
927 
928  /* Set values to zero */
929  for (i = 0; i < 4; i++) e[i] = 0.0;
930 
931  /* Compute values */
932  for (j = -hh ; j <= hh ; j++) {
933  for (i = -hw ; i <= hw ; i++) {
934  diff = *imgdiff++;
935  diffgradx = diff * (*gradx++);
936  diffgrady = diff * (*grady++);
937  e[0] += diffgradx * i + diffgrady * j;
938  e[1] += diffgrady * i - diffgradx * j;
939  e[2] += diffgradx;
940  e[3] += diffgrady;
941  }
942  }
943 
944  for (i = 0; i < 4; i++) e[i] *= 0.5;
945 
946 }
947 
948 void Compute8by8GradientMatrix(KLT_TYPE* gradx,
949  KLT_TYPE* grady,
950  KLT_TYPE a11,
951  KLT_TYPE a12,
952  KLT_TYPE a21,
953  KLT_TYPE a22,
954  KLT_TYPE tx,
955  KLT_TYPE ty,
956  KLT_TYPE h31,
957  KLT_TYPE h32,
958  const int hw,
959  const int hh,
960  Matrix<KLT_TYPE>& T) {
961 
962  register int i, j;
963  T.SetZero();
964 
965  for (j = -hh ; j <= hh ; j++) {
966  for (i = -hw ; i <= hw ; i++) {
967  const double hom = h31*i+h32*j+1.0;
968  const KLT_TYPE gx = *gradx++ / hom;
969  const KLT_TYPE gy = *grady++ / hom;
970  Vector<double> dIdH(8);
971  dIdH[0] = gx*i;
972  dIdH[1] = gy*i;
973  dIdH[2] = gx*j;
974  dIdH[3] = gy*j;
975  dIdH[4] = gx;
976  dIdH[5] = gy;
977  dIdH[6] = -i/hom * (gx*(a11*i+a12*j+tx) + gy*(a21*i+a22*j+ty));
978  dIdH[7] = -j/hom * (gx*(a11*i+a12*j+tx) + gy*(a21*i+a22*j+ty));
979 
980  T += dIdH.OuterProduct(dIdH);
981 
982  }
983  }
984 
985  /*
986  // regularisation:
987  double mean = 0;
988  for (j = 0 ; j < 8; j++) {
989  mean += fabs(T[j][j]);
990  }
991  // make sure we get no worse conditioning than 10^4
992  mean *= 0.0001;
993  for (j = 0 ; j < 8 ; j++) {
994  T[j][j] += mean;
995  }
996  */
997 
998 
999 
1000 
1001 
1002 }
1003 
1004 void Compute8by1ErrorVector(KLT_TYPE* imgdiff,
1005  KLT_TYPE* gradx,
1006  KLT_TYPE* grady,
1007  KLT_TYPE a11,
1008  KLT_TYPE a12,
1009  KLT_TYPE a21,
1010  KLT_TYPE a22,
1011  KLT_TYPE tx,
1012  KLT_TYPE ty,
1013  KLT_TYPE h31,
1014  KLT_TYPE h32,
1015  const int hw,
1016  const int hh,
1017  Vector<KLT_TYPE>& e){
1018  register int i, j;
1019  register KLT_TYPE diffgradx, diffgrady;
1020 
1021  /* Set values to zero */
1022  for(i = 0; i < 8; i++) e[i] = 0.0;
1023 
1024  /* Compute values */
1025  for (j = -hh ; j <= hh ; j++) {
1026  for (i = -hw ; i <= hw ; i++) {
1027 
1028  const double hom = h31*i+h32*j+1.0;
1029  diffgradx = *imgdiff * (*gradx++) / hom;
1030  diffgrady = *imgdiff++ * (*grady++) / hom;
1031 
1032  e[0] += diffgradx*i;
1033  e[1] += diffgrady*i;
1034  e[2] += diffgradx*j;
1035  e[3] += diffgrady*j;
1036  e[4] += diffgradx;
1037  e[5] += diffgrady;
1038  e[6] += -diffgradx*(a11*i+a12*j+tx)/hom*i
1039  -diffgrady*(a21*i+a22*j+ty)/hom*i;
1040  e[7] += -diffgradx*(a11*i+a12*j+tx)/hom*j
1041  -diffgrady*(a21*i+a22*j+ty)/hom*j;
1042 
1043  }
1044  }
1045 
1046  // this is not correct, maybe it helps convergence ???
1047  for(i = 0; i < 8; i++) e[i] *= 0.5;
1048 
1049 
1050 }
1051 }
1052 
1053 
1054 /**************** end affine part ********************************/
1055 
1056 
1057 //////////////////////////////////////////////////////////////////////////
1058 // instantiation
1059 //////////////////////////////////////////////////////////////////////////
1060 namespace BIAS{
1061 template class BIASMatcher2D_EXPORT TrackerBaseAffine<unsigned char>;
1062 template class BIASMatcher2D_EXPORT TrackerBaseAffine<float>;
1063 
1064 // fill in instances as required
1065 #ifdef BUILD_IMAGE_INT
1066 template class TrackerBaseAffine<int>;
1067 #endif
1068 #ifdef BUILD_IMAGE_CHAR
1069 template class TrackerBaseAffine<char>;
1070 #endif
1071 #ifdef BUILD_IMAGE_SHORT
1072 template class TrackerBaseAffine<short>;
1073 #endif
1074 #ifdef BUILD_IMAGE_USHORT
1075 template class TrackerBaseAffine<unsigned short>;
1076 #endif
1077 #ifdef BUILD_IMAGE_UINT
1078 template class TrackerBaseAffine<unsigned int>;
1079 #endif
1080 #ifdef BUILD_IMAGE_DOUBLE
1081 #endif
1082 }
void Compute4by4GradientMatrix(KLT_TYPE *gradx, KLT_TYPE *grady, const int hw, const int hh, BIAS::Matrix< KLT_TYPE > &T)
void Compute6by6GradientMatrix(KLT_TYPE *gradx, KLT_TYPE *grady, const int hw, const int hh, BIAS::Matrix< KLT_TYPE > &T)
computes and holds the singular value decomposition of a rectangular (not necessarily quadratic) Matr...
Definition: SVD.hh:92
Subscript num_cols() const
Definition: cmat.h:320
Vector< double > Solve(const Vector< double > &y) const
Definition: SVD.cpp:135
a 3x3 Matrix describing projective transformations between planes
Definition: HMatrix.hh:39
point slid out of image
virtual void EvaluateResult_(KLT_TYPE &mad, KLT_TYPE &msd, Matrix< KLT_TYPE > &cov)
fill in computed residui from Track_
void ComputeIntensityDifferenceAffine_(const KLT_TYPE &x2, const KLT_TYPE &y2, const KLT_TYPE &Axx, const KLT_TYPE &Ayx, const KLT_TYPE &Axy, const KLT_TYPE &Ayy, const KLT_TYPE &h31, const KLT_TYPE &h32, const int halfwidth, const int halfheight, KLT_TYPE *imgdiff)
TrackerBaseAffine(bool only4params=false, bool useHomography=false)
Matrix< T > & newsize(Subscript M, Subscript N)
Definition: cmat.h:269
int TrackAffine_(KLT_TYPE p1[2], KLT_TYPE p2[2], const Matrix2x2< KLT_TYPE > &p2A, KLT_TYPE result[2], Matrix2x2< KLT_TYPE > &resultA, KLT_TYPE &error, int &iter, Matrix< KLT_TYPE > &Cov, KLT_TYPE &mad, KLT_TYPE &msd, bool SimilarityTransformOnly=false, bool ComputeCovariance=true)
track using affine warp
void GetJacobian(const HomgPoint2D &x, Matrix< double > &Jac) const
returns jacobian of H: R^2 –&gt; R^2
Definition: HMatrix.cpp:190
const T * GetData() const
get the data pointer the member function itself is const (before {..}) because it doesn&#39;t change the ...
Definition: Vector2.hh:236
long int ComputeCovariance(long int NumErrors, long int NumParams, double *Jac, int *Permutation, double &SumOfSquaredErrors, Matrix< double > &Cov)
Compute covariance matrix from Levenberg-Marquardt resulting Jacobian matrix J(x) and permutation mat...
Definition: Minpack.cpp:289
success (error &lt; maxerror)
void SetZero()
Sets all values to zero.
Definition: Matrix.hh:856
unsigned int GetRows() const
Definition: Matrix.hh:202
no spatial structure is present
void Compute8by8GradientMatrix(KLT_TYPE *gradx, KLT_TYPE *grady, KLT_TYPE a11, KLT_TYPE a12, KLT_TYPE a21, KLT_TYPE a22, KLT_TYPE tx, KLT_TYPE ty, KLT_TYPE h31, KLT_TYPE h32, const int hw, const int hh, BIAS::Matrix< KLT_TYPE > &T)
int GaussJordanElimination(BIAS::Matrix< KLT_TYPE > &a, BIAS::Vector< KLT_TYPE > &b)
T * GetData()
get the pointer to the data array of the matrix (for faster direct memeory access) ...
Definition: Matrix.hh:185
static int Save(const std::string &filename, const ImageBase &img, const enum TFileFormat FileFormat=FF_auto, const bool sync=BIAS_DEFAULT_SYNC, const int c_jpeg_quality=BIAS_DEFAULT_IMAGE_QUALITY, const bool forceNewID=BIAS_DEFAULT_FORCENEWID, const bool &writeMetaData=true)
Export image as file using extrnal libs.
Definition: ImageIO.cpp:725
static std::string toString(const T thenumber, int numzeros=DEFAULT_LEADING_ZEROS)
Converts number to string with leading zeros.
Definition: FileHandling.hh:98
void Compute4by1ErrorVector(KLT_TYPE *imgdiff, KLT_TYPE *gradx, KLT_TYPE *grady, const int hw, const int hh, BIAS::Vector< KLT_TYPE > &e)
void Compute6by1ErrorVector(KLT_TYPE *imgdiff, KLT_TYPE *gradx, KLT_TYPE *grady, const int hw, const int hh, BIAS::Vector< KLT_TYPE > &e)
Matrix< double > Invert()
returns pseudoinverse of A = U * S * V^T A^+ = V * S^+ * U^T
Definition: SVD.cpp:214
void Compute8by1ErrorVector(KLT_TYPE *imgdiff, KLT_TYPE *gradx, KLT_TYPE *grady, KLT_TYPE a11, KLT_TYPE a12, KLT_TYPE a21, KLT_TYPE a22, KLT_TYPE tx, KLT_TYPE ty, KLT_TYPE h31, KLT_TYPE h32, const int hw, const int hh, BIAS::Vector< KLT_TYPE > &e)
virtual int Track_(Vector2< KLT_TYPE > &p1, Vector2< KLT_TYPE > &p2, Vector2< KLT_TYPE > &result, KLT_TYPE &error, int &iter, const Matrix2x2< KLT_TYPE > &AffinePred, Matrix2x2< KLT_TYPE > &AffineResult)
Interface of all tracking algorithms implemented in derived classes.
double NormFrobenius() const
Return Frobenius norm = sqrt(trace(A^t * A)).
Definition: Matrix.hh:897
void GetGradientWinAffine_(const KLT_TYPE &x, const KLT_TYPE &y, const KLT_TYPE &Axx, const KLT_TYPE &Ayx, const KLT_TYPE &Axy, const KLT_TYPE &Ayy, const KLT_TYPE &h31, const KLT_TYPE &h32, const int halfwidth, const int halfheight, KLT_TYPE *out_gradx, KLT_TYPE *out_grady)
void SetZero()
zeroes the image
Definition: ImageBase.hh:83
class BIASGeometryBase_EXPORT HomgPoint2D
maxiter is reached and error is above maxerror
const StorageType ** GetImageDataArray() const
overloaded GetImageDataArray() from ImageBase
Definition: Image.hh:153