[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

affine_registration.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2005-2006 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_AFFINE_REGISTRATION_HXX
37 #define VIGRA_AFFINE_REGISTRATION_HXX
38 
39 #include "mathutil.hxx"
40 #include "matrix.hxx"
41 #include "linear_solve.hxx"
42 #include "tinyvector.hxx"
43 #include "splineimageview.hxx"
44 #include "imagecontainer.hxx"
45 #include "multi_shape.hxx"
46 #include "affinegeometry.hxx"
47 
48 #include <cmath>
49 
50 namespace vigra {
51 
52 /** \addtogroup Registration Image Registration
53 
54  Transform different image into a common coordinate system.
55 */
56 //@{
57 
58 /********************************************************/
59 /* */
60 /* affineMatrix2DFromCorrespondingPoints */
61 /* */
62 /********************************************************/
63 
64 /** \brief Create homogeneous matrix that maps corresponding points onto each other.
65 
66  For use with \ref affineWarpImage(). When only two corresponding points are given,
67  the matrix will only represent a similarity transform
68  (translation, rotation, and uniform scaling). When only one point pair is given,
69  the result will be a pure translation.
70 */
71 template <class SrcIterator, class DestIterator>
72 linalg::TemporaryMatrix<double>
73 affineMatrix2DFromCorrespondingPoints(SrcIterator s, SrcIterator send, DestIterator d)
74 {
75  int size = send - s;
76 
77  linalg::TemporaryMatrix<double> ret(identityMatrix<double>(3));
78 
79  if(size == 1)
80  {
81  ret(0,2) = (*d)[0] - (*s)[0];
82  ret(1,2) = (*d)[1] - (*s)[1];
83  }
84  else if(size == 2)
85  {
86  Matrix<double> m(4,4), r(4,1), so(4,1);
87 
88  for(int k=0; k<size; ++k, ++s, ++d)
89  {
90  m(2*k,0) = (*s)[0];
91  m(2*k,1) = -(*s)[1];
92  m(2*k,2) = 1.0;
93  m(2*k,3) = 0.0;
94  r(2*k,0) = (*d)[0];
95 
96  m(2*k+1,0) = (*s)[1];
97  m(2*k+1,1) = (*s)[0];
98  m(2*k+1,2) = 0.0;
99  m(2*k+1,3) = 1.0;
100  r(2*k+1,0) = (*d)[1];
101  }
102 
103  if(!linearSolve(m, r, so))
104  vigra_fail("affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
105 
106  ret(0,0) = so(0,0);
107  ret(1,1) = so(0,0);
108  ret(0,1) = -so(1,0);
109  ret(1,0) = so(1,0);
110  ret(0,2) = so(2,0);
111  ret(1,2) = so(3,0);
112  }
113  else if(size >= 3)
114  {
115  Matrix<double> m(3,3), rx(3,1), sx(3,1), ry(3,1), sy(3,1), c(3,1);
116  c(2,0) = 1.0;
117  for(int k=0; k<size; ++k, ++s, ++d)
118  {
119  c(0,0) = (*s)[0];
120  c(1,0) = (*s)[1];
121 
122  m += outer(c);
123  rx += (*d)[0]*c;
124  ry += (*d)[1]*c;
125  }
126 
127  if(!linearSolve(m, rx, sx) || !linearSolve(m, ry, sy))
128  vigra_fail("affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
129 
130  ret(0,0) = sx(0,0);
131  ret(0,1) = sx(1,0);
132  ret(0,2) = sx(2,0);
133  ret(1,0) = sy(0,0);
134  ret(1,1) = sy(1,0);
135  ret(1,2) = sy(2,0);
136  }
137 
138  return ret;
139 }
140 
141  /** \brief Option object for affine registration functions.
142 
143  The template parameter <tt>SPLINEORDER</tt> (default: 2) specifies
144  the order of interpolation for the intensities at non-integer image
145  coordinates.
146  */
147 template <int SPLINEORDER = 2>
148 class AffineMotionEstimationOptions
149 {
150  public:
151  double burt_filter_strength;
152  int highest_level, iterations_per_level;
153  bool use_laplacian_pyramid;
154 
155  AffineMotionEstimationOptions()
156  : burt_filter_strength(0.4),
157  highest_level(4),
158  iterations_per_level(4),
159  use_laplacian_pyramid(false)
160  {}
161 
162  template <int ORDER>
163  AffineMotionEstimationOptions(AffineMotionEstimationOptions<ORDER> const & other)
164  : burt_filter_strength(other.burt_filter_strength),
165  highest_level(other.highest_level),
166  iterations_per_level(other.iterations_per_level),
167  use_laplacian_pyramid(other.use_laplacian_pyramid)
168  {}
169 
170  /** \brief Change the spline order for intensity interpolation.
171 
172  Usage:
173  \code
174  // use linear interpolation
175  AffineMotionEstimationOptions<>().splineOrder<1>();
176  \endcode
177 
178  Default: order = 2 (quadratic interpolation)
179  */
180  template <int NEWORDER>
181  AffineMotionEstimationOptions<NEWORDER> splineOrder() const
182  {
183  return AffineMotionEstimationOptions<NEWORDER>(*this);
184  }
185 
186  /** \brief Define amount of smoothing before subsampling to the next pyramid level.
187 
188  Pyramids are created with the Burt filter:
189  \code
190  [0.25 - center / 2.0, 0.25, center, 0.25, 0.25 - center / 2.0]
191  \endcode
192  \a center must thus be between 0.25 and 0.5, and the smaller it is,
193  the more smoothing is applied.
194 
195  Default: 0.4 (moderate smoothing)
196  */
197  AffineMotionEstimationOptions & burtFilterCenterStrength(double center)
198  {
199  vigra_precondition(0.25 <= center && center <= 0.5,
200  "AffineMotionEstimationOptions::burtFilterCenterStrength(): center must be between 0.25 and 0.5 (inclusive).");
201  burt_filter_strength = center;
202  return *this;
203  }
204 
205  /** \brief Define the highest level of the image pyramid.
206 
207  The original image is at level 0, and each level downsamples
208  the image by 1/2.
209 
210  Default: 4 (16-fold downsampling)
211  */
212  AffineMotionEstimationOptions & highestPyramidLevel(unsigned int level)
213  {
214  highest_level = (int)level;
215  return *this;
216  }
217 
218  /** \brief Number of Gauss-Newton iterations per level.
219 
220  Default: 4
221  */
222  AffineMotionEstimationOptions & iterationsPerLevel(unsigned int iter)
223  {
224  vigra_precondition(0 < iter,
225  "AffineMotionEstimationOptions::iterationsPerLevel(): must do at least one iteration per level.");
226  iterations_per_level = (int)iter;
227  return *this;
228  }
229 
230  /** \brief Base registration on a Gaussian pyramid.
231 
232  Images are registered such that the similarity in intensity is
233  maximized.
234 
235  Default: true
236  */
237  AffineMotionEstimationOptions & useGaussianPyramid(bool f = true)
238  {
239  use_laplacian_pyramid = !f;
240  return *this;
241  }
242 
243  /** \brief Base registration on a Laplacian pyramid.
244 
245  Images are registered such that the similarity in second
246  derivatives (=Laplacian operator results) is maximized.
247 
248  Default: false
249  */
250  AffineMotionEstimationOptions & useLaplacianPyramid(bool f = true)
251  {
252  use_laplacian_pyramid = f;
253  return *this;
254  }
255 };
256 
257 namespace detail {
258 
259 struct TranslationEstimationFunctor
260 {
261  template <class SplineImage, class Image>
262  void operator()(SplineImage const & src, Image const & dest, Matrix<double> & matrix) const
263  {
264  int w = dest.width();
265  int h = dest.height();
266 
267  Matrix<double> grad(2,1), m(2,2), r(2,1), s(2,1);
268  double dx = matrix(0,0), dy = matrix(1,0);
269 
270  for(int y = 0; y < h; ++y)
271  {
272  double sx = matrix(0,1)*y + matrix(0,2);
273  double sy = matrix(1,1)*y + matrix(1,2);
274  for(int x = 0; x < w; ++x, sx += dx, sy += dy)
275  {
276  if(!src.isInside(sx, sy))
277  continue;
278 
279  grad(0,0) = src.dx(sx, sy);
280  grad(1,0) = src.dy(sx, sy);
281  double diff = dest(x, y) - src(sx, sy);
282 
283  m += outer(grad);
284  r -= diff*grad;
285  }
286  }
287 
288  linearSolve(m, r, s);
289 
290  matrix(0,2) -= s(0,0);
291  matrix(1,2) -= s(1,0);
292  }
293 };
294 
295 struct SimilarityTransformEstimationFunctor
296 {
297  template <class SplineImage, class Image>
298  void operator()(SplineImage const & src, Image const & dest, Matrix<double> & matrix) const
299  {
300  int w = dest.width();
301  int h = dest.height();
302 
303  Matrix<double> grad(2,1), coord(4, 2), c(4, 1), m(4, 4), r(4,1), s(4,1);
304  coord(0,0) = 1.0;
305  coord(1,1) = 1.0;
306  double dx = matrix(0,0), dy = matrix(1,0);
307 
308  for(int y = 0; y < h; ++y)
309  {
310  double sx = matrix(0,1)*y + matrix(0,2);
311  double sy = matrix(1,1)*y + matrix(1,2);
312  for(int x = 0; x < w; ++x, sx += dx, sy += dy)
313  {
314  if(!src.isInside(sx, sy))
315  continue;
316 
317  grad(0,0) = src.dx(sx, sy);
318  grad(1,0) = src.dy(sx, sy);
319  coord(2,0) = (double)x;
320  coord(3,1) = (double)x;
321  coord(3,0) = -(double)y;
322  coord(2,1) = (double)y;
323  double diff = dest(x, y) - src(sx, sy);
324 
325  c = coord * grad;
326  m += outer(c);
327  r -= diff*c;
328  }
329  }
330 
331  linearSolve(m, r, s);
332 
333  matrix(0,2) -= s(0,0);
334  matrix(1,2) -= s(1,0);
335  matrix(0,0) -= s(2,0);
336  matrix(1,1) -= s(2,0);
337  matrix(0,1) += s(3,0);
338  matrix(1,0) -= s(3,0);
339  }
340 };
341 
342 struct AffineTransformEstimationFunctor
343 {
344  template <class SplineImage, class Image>
345  void operator()(SplineImage const & src, Image const & dest, Matrix<double> & matrix) const
346  {
347  int w = dest.width();
348  int h = dest.height();
349 
350  Matrix<double> grad(2,1), coord(6, 2), c(6, 1), m(6,6), r(6,1), s(6,1);
351  coord(0,0) = 1.0;
352  coord(1,1) = 1.0;
353  double dx = matrix(0,0), dy = matrix(1,0);
354 
355  for(int y = 0; y < h; ++y)
356  {
357  double sx = matrix(0,1)*y + matrix(0,2);
358  double sy = matrix(1,1)*y + matrix(1,2);
359  for(int x = 0; x < w; ++x, sx += dx, sy += dy)
360  {
361  if(!src.isInside(sx, sy))
362  continue;
363 
364  grad(0,0) = src.dx(sx, sy);
365  grad(1,0) = src.dy(sx, sy);
366  coord(2,0) = (double)x;
367  coord(4,1) = (double)x;
368  coord(3,0) = (double)y;
369  coord(5,1) = (double)y;
370  double diff = dest(x, y) - src(sx, sy);
371 
372  c = coord * grad;
373  m += outer(c);
374  r -= diff*c;
375  }
376  }
377 
378  linearSolve(m, r, s);
379 
380  matrix(0,2) -= s(0,0);
381  matrix(1,2) -= s(1,0);
382  matrix(0,0) -= s(2,0);
383  matrix(0,1) -= s(3,0);
384  matrix(1,0) -= s(4,0);
385  matrix(1,1) -= s(5,0);
386  }
387 };
388 
389 template <class SrcIterator, class SrcAccessor,
390  class DestIterator, class DestAccessor,
391  int SPLINEORDER, class Functor>
392 void
393 estimateAffineMotionImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
394  DestIterator dul, DestIterator dlr, DestAccessor dest,
395  Matrix<double> & affineMatrix,
396  AffineMotionEstimationOptions<SPLINEORDER> const & options,
397  Functor motionModel)
398 {
399  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote STmpType;
400  typedef BasicImage<STmpType> STmpImage;
401  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote DTmpType;
402  typedef BasicImage<DTmpType> DTmpImage;
403 
404  int toplevel = options.highest_level;
405  ImagePyramid<STmpImage> srcPyramid(0, toplevel, sul, slr, src);
406  ImagePyramid<DTmpImage> destPyramid(0, toplevel, dul, dlr, dest);
407 
408  if(options.use_laplacian_pyramid)
409  {
410  pyramidReduceBurtLaplacian(srcPyramid, 0, toplevel, options.burt_filter_strength);
411  pyramidReduceBurtLaplacian(destPyramid, 0, toplevel, options.burt_filter_strength);
412  }
413  else
414  {
415  pyramidReduceBurtFilter(srcPyramid, 0, toplevel, options.burt_filter_strength);
416  pyramidReduceBurtFilter(destPyramid, 0, toplevel, options.burt_filter_strength);
417  }
418 
419  Matrix<double> currentMatrix(affineMatrix(2,2) == 0.0
420  ? identityMatrix<double>(3)
421  : affineMatrix);
422  currentMatrix(0,2) /= std::pow(2.0, toplevel);
423  currentMatrix(1,2) /= std::pow(2.0, toplevel);
424 
425  for(int level = toplevel; level >= 0; --level)
426  {
427  SplineImageView<SPLINEORDER, STmpType> sp(srcImageRange(srcPyramid[level]));
428 
429  for(int iter = 0; iter < options.iterations_per_level; ++iter)
430  {
431  motionModel(sp, destPyramid[level], currentMatrix);
432  }
433 
434  if(level > 0)
435  {
436  currentMatrix(0,2) *= 2.0;
437  currentMatrix(1,2) *= 2.0;
438  }
439  }
440 
441  affineMatrix = currentMatrix;
442 }
443 
444 } // namespace detail
445 
446 /********************************************************/
447 /* */
448 /* estimateTranslation */
449 /* */
450 /********************************************************/
451 
452 /** \brief Estimate the optical flow between two images according to a translation model.
453 
454  This function applies the same algorithm as \ref estimateAffineTransform()
455  with the additional constraint that the motion model must be a translation
456  rather than affine.
457 
458  <b> Declarations:</b>
459 
460  <b>\#include</b> <vigra/affine_registration.hxx><br>
461  Namespace: vigra
462 
463  pass 2D array views:
464  \code
465  namespace vigra {
466  template <class T1, class S1,
467  class T2, class S2,
468  int SPLINEORDER>
469  void
470  estimateTranslation(MultiArrayView<2, T1, S1> const & src,
471  MultiArrayView<2, T2, S2> dest,
472  Matrix<double> & affineMatrix,
473  AffineMotionEstimationOptions<SPLINEORDER> const & options = AffineMotionEstimationOptions<SPLINEORDER>());
474  }
475  \endcode
476 
477  \deprecatedAPI{estimateTranslation}
478  pass \ref ImageIterators and \ref DataAccessors :
479  \code
480  namespace vigra {
481  template <class SrcIterator, class SrcAccessor,
482  class DestIterator, class DestAccessor,
483  int SPLINEORDER = 2>
484  void
485  estimateTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
486  DestIterator dul, DestIterator dlr, DestAccessor dest,
487  Matrix<double> & affineMatrix,
488  AffineMotionEstimationOptions<SPLINEORDER> const & options =
489  AffineMotionEstimationOptions<>())
490  }
491  \endcode
492  use argument objects in conjunction with \ref ArgumentObjectFactories :
493  \code
494  namespace vigra {
495  template <class SrcIterator, class SrcAccessor,
496  class DestIterator, class DestAccessor,
497  int SPLINEORDER = 2>
498  void
499  estimateTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
500  triple<DestIterator, DestIterator, DestAccessor> dest,
501  Matrix<double> & affineMatrix,
502  AffineMotionEstimationOptions<SPLINEORDER> const & options =
503  AffineMotionEstimationOptions<>())
504  }
505  \endcode
506  \deprecatedEnd
507 */
508 doxygen_overloaded_function(template <...> void estimateTranslation)
509 
510 template <class SrcIterator, class SrcAccessor,
511  class DestIterator, class DestAccessor,
512  int SPLINEORDER>
513 inline void
514 estimateTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
515  DestIterator dul, DestIterator dlr, DestAccessor dest,
516  Matrix<double> & affineMatrix,
517  AffineMotionEstimationOptions<SPLINEORDER> const & options)
518 {
519  detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
520  options, detail::TranslationEstimationFunctor());
521 }
522 
523 template <class SrcIterator, class SrcAccessor,
524  class DestIterator, class DestAccessor>
525 inline void
526 estimateTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
527  DestIterator dul, DestIterator dlr, DestAccessor dest,
528  Matrix<double> & affineMatrix)
529 {
530  estimateTranslation(sul, slr, src, dul, dlr, dest,
531  affineMatrix, AffineMotionEstimationOptions<>());
532 }
533 
534 template <class SrcIterator, class SrcAccessor,
535  class DestIterator, class DestAccessor,
536  int SPLINEORDER>
537 inline void
538 estimateTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
539  triple<DestIterator, DestIterator, DestAccessor> dest,
540  Matrix<double> & affineMatrix,
541  AffineMotionEstimationOptions<SPLINEORDER> const & options)
542 {
543  estimateTranslation(src.first, src.second, src.third, dest.first, dest.second, dest.third,
544  affineMatrix, options);
545 }
546 
547 template <class SrcIterator, class SrcAccessor,
548  class DestIterator, class DestAccessor>
549 inline void
550 estimateTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
551  triple<DestIterator, DestIterator, DestAccessor> dest,
552  Matrix<double> & affineMatrix)
553 {
554  estimateTranslation(src.first, src.second, src.third, dest.first, dest.second, dest.third,
555  affineMatrix, AffineMotionEstimationOptions<>());
556 }
557 
558 template <class T1, class S1,
559  class T2, class S2,
560  int SPLINEORDER>
561 inline void
562 estimateTranslation(MultiArrayView<2, T1, S1> const & src,
563  MultiArrayView<2, T2, S2> dest,
564  Matrix<double> & affineMatrix,
565  AffineMotionEstimationOptions<SPLINEORDER> const & options)
566 {
567  estimateTranslation(srcImageRange(src), destImageRange(dest),
568  affineMatrix, options);
569 }
570 
571 template <class T1, class S1,
572  class T2, class S2>
573 inline void
574 estimateTranslation(MultiArrayView<2, T1, S1> const & src,
575  MultiArrayView<2, T2, S2> dest,
576  Matrix<double> & affineMatrix)
577 {
578  estimateTranslation(srcImageRange(src), destImageRange(dest),
579  affineMatrix, AffineMotionEstimationOptions<>());
580 }
581 
582 /********************************************************/
583 /* */
584 /* estimateSimilarityTransform */
585 /* */
586 /********************************************************/
587 
588 /** \brief Estimate the optical flow between two images according to a similarity transform model
589  (e.g. translation, rotation, and uniform scaling).
590 
591  This function applies the same algorithm as \ref estimateAffineTransform()
592  with the additional constraint that the motion model must be a similarity
593  transform rather than affine.
594 
595  <b> Declarations:</b>
596 
597  <b>\#include</b> <vigra/affine_registration.hxx><br>
598  Namespace: vigra
599 
600  pass 2D array views:
601  \code
602  namespace vigra {
603  template <class T1, class S1,
604  class T2, class S2,
605  int SPLINEORDER>
606  void
607  estimateSimilarityTransform(MultiArrayView<2, T1, S1> const & src,
608  MultiArrayView<2, T2, S2> dest,
609  Matrix<double> & affineMatrix,
610  AffineMotionEstimationOptions<SPLINEORDER> const & options =
611  AffineMotionEstimationOptions<SPLINEORDER>());
612  }
613  \endcode
614 
615  \deprecatedAPI{estimateSimilarityTransform}
616  pass \ref ImageIterators and \ref DataAccessors :
617  \code
618  namespace vigra {
619  template <class SrcIterator, class SrcAccessor,
620  class DestIterator, class DestAccessor,
621  int SPLINEORDER = 2>
622  void
623  estimateSimilarityTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
624  DestIterator dul, DestIterator dlr, DestAccessor dest,
625  Matrix<double> & affineMatrix,
626  AffineMotionEstimationOptions<SPLINEORDER> const & options =
627  AffineMotionEstimationOptions<>())
628  }
629  \endcode
630  use argument objects in conjunction with \ref ArgumentObjectFactories :
631  \code
632  namespace vigra {
633  template <class SrcIterator, class SrcAccessor,
634  class DestIterator, class DestAccessor,
635  int SPLINEORDER = 2>
636  void
637  estimateSimilarityTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
638  triple<DestIterator, DestIterator, DestAccessor> dest,
639  Matrix<double> & affineMatrix,
640  AffineMotionEstimationOptions<SPLINEORDER> const & options =
641  AffineMotionEstimationOptions<>())
642  }
643  \endcode
644  \deprecatedEnd
645 */
646 doxygen_overloaded_function(template <...> void estimateSimilarityTransform)
647 
648 template <class SrcIterator, class SrcAccessor,
649  class DestIterator, class DestAccessor,
650  int SPLINEORDER>
651 inline void
652 estimateSimilarityTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
653  DestIterator dul, DestIterator dlr, DestAccessor dest,
654  Matrix<double> & affineMatrix,
655  AffineMotionEstimationOptions<SPLINEORDER> const & options)
656 {
657  detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
658  options, detail::SimilarityTransformEstimationFunctor());
659 }
660 
661 template <class SrcIterator, class SrcAccessor,
662  class DestIterator, class DestAccessor>
663 inline void
664 estimateSimilarityTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
665  DestIterator dul, DestIterator dlr, DestAccessor dest,
666  Matrix<double> & affineMatrix)
667 {
668  estimateSimilarityTransform(sul, slr, src, dul, dlr, dest,
669  affineMatrix, AffineMotionEstimationOptions<>());
670 }
671 
672 template <class SrcIterator, class SrcAccessor,
673  class DestIterator, class DestAccessor,
674  int SPLINEORDER>
675 inline void
676 estimateSimilarityTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
677  triple<DestIterator, DestIterator, DestAccessor> dest,
678  Matrix<double> & affineMatrix,
679  AffineMotionEstimationOptions<SPLINEORDER> const & options)
680 {
681  estimateSimilarityTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
682  affineMatrix, options);
683 }
684 
685 template <class SrcIterator, class SrcAccessor,
686  class DestIterator, class DestAccessor>
687 inline void
688 estimateSimilarityTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
689  triple<DestIterator, DestIterator, DestAccessor> dest,
690  Matrix<double> & affineMatrix)
691 {
692  estimateSimilarityTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
693  affineMatrix, AffineMotionEstimationOptions<>());
694 }
695 
696 template <class T1, class S1,
697  class T2, class S2,
698  int SPLINEORDER>
699 inline void
700 estimateSimilarityTransform(MultiArrayView<2, T1, S1> const & src,
701  MultiArrayView<2, T2, S2> dest,
702  Matrix<double> & affineMatrix,
703  AffineMotionEstimationOptions<SPLINEORDER> const & options)
704 {
705  estimateSimilarityTransform(srcImageRange(src), destImageRange(dest),
706  affineMatrix, options);
707 }
708 
709 template <class T1, class S1,
710  class T2, class S2>
711 inline void
712 estimateSimilarityTransform(MultiArrayView<2, T1, S1> const & src,
713  MultiArrayView<2, T2, S2> dest,
714  Matrix<double> & affineMatrix)
715 {
716  estimateSimilarityTransform(srcImageRange(src), destImageRange(dest),
717  affineMatrix, AffineMotionEstimationOptions<>());
718 }
719 
720 /********************************************************/
721 /* */
722 /* estimateAffineTransform */
723 /* */
724 /********************************************************/
725 
726 /** \brief Estimate the optical flow between two images according to an affine transform model
727  (e.g. translation, rotation, non-uniform scaling, and shearing).
728 
729  This function implements the algorithm described in
730 
731  J.R. Bergen, P. Anandan, K.J. Hanna, R. Hingorani: <i>"Hierarchical model-based motion estimation"</i>, ECCV 1992
732 
733  Specifically, it minimizes the squared loss between the images I at two consecutive time
734  points t-1 and t:
735  \f[ \min_{\theta} \sum_{\mathbf{x}} \left(I(\mathbf{x}, t) - I(\mathbf{x} - \mathbf{u}_{\theta}(\mathbf{x}), t-1)\right)^2
736  \f]
737  where \f$\mathbf{x}\f$ are the pixel coordinates and \f$\mathbf{u}_{\theta}(\mathbf{x})\f$
738  is an affine motion model parameterized by \f$\theta\f$. Since the objective is
739  non-linear, it is linearized by first-order Taylor expansion w.r.t. \f$\theta\f$,
740  and a local optimum is determined iteratively by the Gauss-Newton method. To handle
741  larger displacements, the algorithm employs a coarse-to-fine strategy, where the
742  motion is first estimated on downsampled versions of the images and then refined at
743  consecutively higher resolutions.
744 
745  The algorithm's parameters can be controlled by the option object
746  \ref vigra::AffineMotionEstimationOptions. In particular, one can determine if
747  <ul>
748  <li> I in the objective refers to the original images (default) or their second
749  derivatives (<tt>options.useLaplacianPyramid()</tt> -- makes motion
750  estimation invariant against additive intensity offsets);</li>
751  <li> the highest pyramid level to be used
752  (<tt>options.highestPyramidLevel(h)</tt> -- images are downsampled to 2<sup>-h</sup> times their original size, default: h=4);</li>
753  <li> the number of Gauss-Newton iterations per resolution level
754  (<tt>options.iterationsPerLevel(i)</tt>, default: i=4);</li>
755  <li> the interpolation order to compute subpixel intensities \f$I(\mathbf{x} - \mathbf{u}_{\theta}(\mathbf{x}), t-1)\f$ (default: 2)
756  </ul>
757 
758  The resulting affine model is stored in parameter <tt>affineMatrix</tt>, which
759  can be used by \ref affineWarpImage() to apply the transformation to time frame t-1.
760  See documentation there for the precise meaning of the matrix elements.
761 
762  <b> Declarations:</b>
763 
764  <b>\#include</b> <vigra/affine_registration.hxx><br>
765  Namespace: vigra
766 
767  pass 2D array views:
768  \code
769  namespace vigra {
770  template <class T1, class S1,
771  class T2, class S2,
772  int SPLINEORDER>
773  void
774  estimateAffineTransform(MultiArrayView<2, T1, S1> const & src,
775  MultiArrayView<2, T2, S2> dest,
776  Matrix<double> & affineMatrix,
777  AffineMotionEstimationOptions<SPLINEORDER> const & options =
778  AffineMotionEstimationOptions<SPLINEORDER>());
779  }
780  \endcode
781 
782  \deprecatedAPI{estimateAffineTransform}
783  pass \ref ImageIterators and \ref DataAccessors :
784  \code
785  namespace vigra {
786  template <class SrcIterator, class SrcAccessor,
787  class DestIterator, class DestAccessor,
788  int SPLINEORDER = 2>
789  void
790  estimateAffineTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
791  DestIterator dul, DestIterator dlr, DestAccessor dest,
792  Matrix<double> & affineMatrix,
793  AffineMotionEstimationOptions<SPLINEORDER> const & options =
794  AffineMotionEstimationOptions<>())
795  }
796  \endcode
797  use argument objects in conjunction with \ref ArgumentObjectFactories :
798  \code
799  namespace vigra {
800  template <class SrcIterator, class SrcAccessor,
801  class DestIterator, class DestAccessor,
802  int SPLINEORDER = 2>
803  void
804  estimateAffineTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
805  triple<DestIterator, DestIterator, DestAccessor> dest,
806  Matrix<double> & affineMatrix,
807  AffineMotionEstimationOptions<SPLINEORDER> const & options =
808  AffineMotionEstimationOptions<>())
809  }
810  \endcode
811  \deprecatedEnd
812 */
813 doxygen_overloaded_function(template <...> void estimateAffineTransform)
814 
815 template <class SrcIterator, class SrcAccessor,
816  class DestIterator, class DestAccessor,
817  int SPLINEORDER>
818 inline void
819 estimateAffineTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
820  DestIterator dul, DestIterator dlr, DestAccessor dest,
821  Matrix<double> & affineMatrix,
822  AffineMotionEstimationOptions<SPLINEORDER> const & options)
823 {
824  detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
825  options, detail::AffineTransformEstimationFunctor());
826 }
827 
828 template <class SrcIterator, class SrcAccessor,
829  class DestIterator, class DestAccessor>
830 inline void
831 estimateAffineTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
832  DestIterator dul, DestIterator dlr, DestAccessor dest,
833  Matrix<double> & affineMatrix)
834 {
835  estimateAffineTransform(sul, slr, src, dul, dlr, dest,
836  affineMatrix, AffineMotionEstimationOptions<>());
837 }
838 
839 template <class SrcIterator, class SrcAccessor,
840  class DestIterator, class DestAccessor,
841  int SPLINEORDER>
842 inline void
843 estimateAffineTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
844  triple<DestIterator, DestIterator, DestAccessor> dest,
845  Matrix<double> & affineMatrix,
846  AffineMotionEstimationOptions<SPLINEORDER> const & options)
847 {
848  estimateAffineTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
849  affineMatrix, options);
850 }
851 
852 template <class SrcIterator, class SrcAccessor,
853  class DestIterator, class DestAccessor>
854 inline void
855 estimateAffineTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
856  triple<DestIterator, DestIterator, DestAccessor> dest,
857  Matrix<double> & affineMatrix)
858 {
859  estimateAffineTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
860  affineMatrix, AffineMotionEstimationOptions<>());
861 }
862 
863 template <class T1, class S1,
864  class T2, class S2,
865  int SPLINEORDER>
866 inline void
867 estimateAffineTransform(MultiArrayView<2, T1, S1> const & src,
868  MultiArrayView<2, T2, S2> dest,
869  Matrix<double> & affineMatrix,
870  AffineMotionEstimationOptions<SPLINEORDER> const & options)
871 {
872  vigra_precondition(src.shape() == dest.shape(),
873  "estimateAffineTransform(): shape mismatch between input and output.");
874  estimateAffineTransform(srcImageRange(src), destImageRange(dest),
875  affineMatrix, options);
876 }
877 
878 template <class T1, class S1,
879  class T2, class S2>
880 inline void
881 estimateAffineTransform(MultiArrayView<2, T1, S1> const & src,
882  MultiArrayView<2, T2, S2> dest,
883  Matrix<double> & affineMatrix)
884 {
885  vigra_precondition(src.shape() == dest.shape(),
886  "estimateAffineTransform(): shape mismatch between input and output.");
887  estimateAffineTransform(srcImageRange(src), destImageRange(dest),
888  affineMatrix, AffineMotionEstimationOptions<>());
889 }
890 
891 //@}
892 
893 } // namespace vigra
894 
895 
896 #endif /* VIGRA_AFFINE_REGISTRATION_HXX */
vigra::AffineMotionEstimationOptions::highestPyramidLevel
AffineMotionEstimationOptions & highestPyramidLevel(unsigned int level)
Define the highest level of the image pyramid.
Definition: affine_registration.hxx:212
vigra::AffineMotionEstimationOptions::useGaussianPyramid
AffineMotionEstimationOptions & useGaussianPyramid(bool f=true)
Base registration on a Gaussian pyramid.
Definition: affine_registration.hxx:237
vigra::AffineMotionEstimationOptions::splineOrder
AffineMotionEstimationOptions< NEWORDER > splineOrder() const
Change the spline order for intensity interpolation.
Definition: affine_registration.hxx:181
vigra::AffineMotionEstimationOptions
Option object for affine registration functions.
Definition: affine_registration.hxx:148
vigra::AffineMotionEstimationOptions::iterationsPerLevel
AffineMotionEstimationOptions & iterationsPerLevel(unsigned int iter)
Number of Gauss-Newton iterations per level.
Definition: affine_registration.hxx:222
vigra::AffineMotionEstimationOptions::burtFilterCenterStrength
AffineMotionEstimationOptions & burtFilterCenterStrength(double center)
Define amount of smoothing before subsampling to the next pyramid level.
Definition: affine_registration.hxx:197
vigra::estimateSimilarityTransform
void estimateSimilarityTransform(...)
Estimate the optical flow between two images according to a similarity transform model (e...
vigra::estimateAffineTransform
void estimateAffineTransform(...)
Estimate the optical flow between two images according to an affine transform model (e...
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::estimateTranslation
void estimateTranslation(...)
Estimate the optical flow between two images according to a translation model.
vigra::pyramidReduceBurtFilter
void pyramidReduceBurtFilter(...)
Two-fold down-sampling for image pyramid construction.
vigra::linalg::outer
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1459
vigra::AffineMotionEstimationOptions::useLaplacianPyramid
AffineMotionEstimationOptions & useLaplacianPyramid(bool f=true)
Base registration on a Laplacian pyramid.
Definition: affine_registration.hxx:250
vigra::pyramidReduceBurtLaplacian
void pyramidReduceBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Create a Laplacian pyramid.
Definition: resampling_convolution.hxx:1197
vigra::linalg::linearSolve
bool linearSolve(...)
vigra::affineMatrix2DFromCorrespondingPoints
linalg::TemporaryMatrix< double > affineMatrix2DFromCorrespondingPoints(SrcIterator s, SrcIterator send, DestIterator d)
Create homogeneous matrix that maps corresponding points onto each other.
Definition: affine_registration.hxx:73

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)