36 #ifndef VIGRA_AFFINE_REGISTRATION_HXX
37 #define VIGRA_AFFINE_REGISTRATION_HXX
39 #include "mathutil.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"
71 template <
class SrcIterator,
class DestIterator>
72 linalg::TemporaryMatrix<double>
77 linalg::TemporaryMatrix<double> ret(identityMatrix<double>(3));
81 ret(0,2) = (*d)[0] - (*s)[0];
82 ret(1,2) = (*d)[1] - (*s)[1];
88 for(
int k=0; k<size; ++k, ++s, ++d)
100 r(2*k+1,0) = (*d)[1];
103 if(!linearSolve(m, r, so))
104 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
115 Matrix<double> m(3,3), rx(3,1), sx(3,1), ry(3,1), sy(3,1), c(3,1);
117 for(
int k=0; k<size; ++k, ++s, ++d)
127 if(!linearSolve(m, rx, sx) || !linearSolve(m, ry, sy))
128 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
147 template <
int SPLINEORDER = 2>
151 double burt_filter_strength;
152 int highest_level, iterations_per_level;
153 bool use_laplacian_pyramid;
156 : burt_filter_strength(0.4),
158 iterations_per_level(4),
159 use_laplacian_pyramid(
false)
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)
180 template <
int NEWORDER>
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;
214 highest_level = (int)level;
224 vigra_precondition(0 < iter,
225 "AffineMotionEstimationOptions::iterationsPerLevel(): must do at least one iteration per level.");
226 iterations_per_level = (int)iter;
239 use_laplacian_pyramid = !f;
252 use_laplacian_pyramid = f;
259 struct TranslationEstimationFunctor
261 template <
class SplineImage,
class Image>
262 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const
264 int w = dest.width();
265 int h = dest.height();
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);
270 for(
int y = 0; y < h; ++y)
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)
276 if(!src.isInside(sx, sy))
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);
290 matrix(0,2) -= s(0,0);
291 matrix(1,2) -= s(1,0);
295 struct SimilarityTransformEstimationFunctor
297 template <
class SplineImage,
class Image>
298 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const
300 int w = dest.width();
301 int h = dest.height();
303 Matrix<double> grad(2,1), coord(4, 2), c(4, 1), m(4, 4), r(4,1), s(4,1);
306 double dx = matrix(0,0), dy = matrix(1,0);
308 for(
int y = 0; y < h; ++y)
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)
314 if(!src.isInside(sx, sy))
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);
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);
342 struct AffineTransformEstimationFunctor
344 template <
class SplineImage,
class Image>
345 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const
347 int w = dest.width();
348 int h = dest.height();
350 Matrix<double> grad(2,1), coord(6, 2), c(6, 1), m(6,6), r(6,1), s(6,1);
353 double dx = matrix(0,0), dy = matrix(1,0);
355 for(
int y = 0; y < h; ++y)
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)
361 if(!src.isInside(sx, sy))
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);
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);
389 template <
class SrcIterator,
class SrcAccessor,
390 class DestIterator,
class DestAccessor,
391 int SPLINEORDER,
class Functor>
393 estimateAffineMotionImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
394 DestIterator dul, DestIterator dlr, DestAccessor dest,
395 Matrix<double> & affineMatrix,
396 AffineMotionEstimationOptions<SPLINEORDER>
const & options,
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;
404 int toplevel = options.highest_level;
405 ImagePyramid<STmpImage> srcPyramid(0, toplevel, sul, slr, src);
406 ImagePyramid<DTmpImage> destPyramid(0, toplevel, dul, dlr, dest);
408 if(options.use_laplacian_pyramid)
419 Matrix<double> currentMatrix(affineMatrix(2,2) == 0.0
420 ? identityMatrix<double>(3)
422 currentMatrix(0,2) /= std::pow(2.0, toplevel);
423 currentMatrix(1,2) /= std::pow(2.0, toplevel);
425 for(
int level = toplevel; level >= 0; --level)
427 SplineImageView<SPLINEORDER, STmpType> sp(srcImageRange(srcPyramid[level]));
429 for(
int iter = 0; iter < options.iterations_per_level; ++iter)
431 motionModel(sp, destPyramid[level], currentMatrix);
436 currentMatrix(0,2) *= 2.0;
437 currentMatrix(1,2) *= 2.0;
441 affineMatrix = currentMatrix;
510 template <
class SrcIterator,
class SrcAccessor,
511 class DestIterator,
class DestAccessor,
515 DestIterator dul, DestIterator dlr, DestAccessor dest,
519 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
520 options, detail::TranslationEstimationFunctor());
523 template <
class SrcIterator,
class SrcAccessor,
524 class DestIterator,
class DestAccessor>
527 DestIterator dul, DestIterator dlr, DestAccessor dest,
528 Matrix<double> & affineMatrix)
531 affineMatrix, AffineMotionEstimationOptions<>());
534 template <
class SrcIterator,
class SrcAccessor,
535 class DestIterator,
class DestAccessor,
539 triple<DestIterator, DestIterator, DestAccessor> dest,
540 Matrix<double> & affineMatrix,
541 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
544 affineMatrix, options);
547 template <
class SrcIterator,
class SrcAccessor,
548 class DestIterator,
class DestAccessor>
551 triple<DestIterator, DestIterator, DestAccessor> dest,
552 Matrix<double> & affineMatrix)
555 affineMatrix, AffineMotionEstimationOptions<>());
558 template <
class T1,
class S1,
563 MultiArrayView<2, T2, S2> dest,
564 Matrix<double> & affineMatrix,
565 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
568 affineMatrix, options);
571 template <
class T1,
class S1,
575 MultiArrayView<2, T2, S2> dest,
576 Matrix<double> & affineMatrix)
579 affineMatrix, AffineMotionEstimationOptions<>());
648 template <
class SrcIterator,
class SrcAccessor,
649 class DestIterator,
class DestAccessor,
653 DestIterator dul, DestIterator dlr, DestAccessor dest,
657 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
658 options, detail::SimilarityTransformEstimationFunctor());
661 template <
class SrcIterator,
class SrcAccessor,
662 class DestIterator,
class DestAccessor>
665 DestIterator dul, DestIterator dlr, DestAccessor dest,
666 Matrix<double> & affineMatrix)
669 affineMatrix, AffineMotionEstimationOptions<>());
672 template <
class SrcIterator,
class SrcAccessor,
673 class DestIterator,
class DestAccessor,
677 triple<DestIterator, DestIterator, DestAccessor> dest,
678 Matrix<double> & affineMatrix,
679 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
682 affineMatrix, options);
685 template <
class SrcIterator,
class SrcAccessor,
686 class DestIterator,
class DestAccessor>
689 triple<DestIterator, DestIterator, DestAccessor> dest,
690 Matrix<double> & affineMatrix)
693 affineMatrix, AffineMotionEstimationOptions<>());
696 template <
class T1,
class S1,
701 MultiArrayView<2, T2, S2> dest,
702 Matrix<double> & affineMatrix,
703 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
706 affineMatrix, options);
709 template <
class T1,
class S1,
713 MultiArrayView<2, T2, S2> dest,
714 Matrix<double> & affineMatrix)
717 affineMatrix, AffineMotionEstimationOptions<>());
815 template <
class SrcIterator,
class SrcAccessor,
816 class DestIterator,
class DestAccessor,
820 DestIterator dul, DestIterator dlr, DestAccessor dest,
824 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
825 options, detail::AffineTransformEstimationFunctor());
828 template <
class SrcIterator,
class SrcAccessor,
829 class DestIterator,
class DestAccessor>
832 DestIterator dul, DestIterator dlr, DestAccessor dest,
833 Matrix<double> & affineMatrix)
836 affineMatrix, AffineMotionEstimationOptions<>());
839 template <
class SrcIterator,
class SrcAccessor,
840 class DestIterator,
class DestAccessor,
844 triple<DestIterator, DestIterator, DestAccessor> dest,
845 Matrix<double> & affineMatrix,
846 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
849 affineMatrix, options);
852 template <
class SrcIterator,
class SrcAccessor,
853 class DestIterator,
class DestAccessor>
856 triple<DestIterator, DestIterator, DestAccessor> dest,
857 Matrix<double> & affineMatrix)
860 affineMatrix, AffineMotionEstimationOptions<>());
863 template <
class T1,
class S1,
868 MultiArrayView<2, T2, S2> dest,
869 Matrix<double> & affineMatrix,
870 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
872 vigra_precondition(src.shape() == dest.shape(),
873 "estimateAffineTransform(): shape mismatch between input and output.");
875 affineMatrix, options);
878 template <
class T1,
class S1,
882 MultiArrayView<2, T2, S2> dest,
883 Matrix<double> & affineMatrix)
885 vigra_precondition(src.shape() == dest.shape(),
886 "estimateAffineTransform(): shape mismatch between input and output.");
888 affineMatrix, AffineMotionEstimationOptions<>());