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" 68 template <
class SrcIterator,
class DestIterator>
69 linalg::TemporaryMatrix<double>
74 linalg::TemporaryMatrix<double> ret(identityMatrix<double>(3));
78 ret(0,2) = (*d)[0] - (*s)[0];
79 ret(1,2) = (*d)[1] - (*s)[1];
83 Matrix<double> m(4,4), r(4,1), so(4,1);
85 for(
int k=0; k<size; ++k, ++s, ++d)
101 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
112 Matrix<double> m(3,3), rx(3,1), sx(3,1), ry(3,1), sy(3,1), c(3,1);
114 for(
int k=0; k<size; ++k, ++s, ++d)
125 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
138 template <
int SPLINEORDER = 2>
139 class AffineMotionEstimationOptions
142 double burt_filter_strength;
143 int highest_level, iterations_per_level;
144 bool use_laplacian_pyramid;
146 AffineMotionEstimationOptions()
147 : burt_filter_strength(0.4),
149 iterations_per_level(4),
150 use_laplacian_pyramid(
false)
154 AffineMotionEstimationOptions(AffineMotionEstimationOptions<ORDER>
const & other)
155 : burt_filter_strength(other.burt_filter_strength),
156 highest_level(other.highest_level),
157 iterations_per_level(other.iterations_per_level),
158 use_laplacian_pyramid(other.use_laplacian_pyramid)
161 template <
int NEWORDER>
162 AffineMotionEstimationOptions<NEWORDER> splineOrder()
const 164 return AffineMotionEstimationOptions<NEWORDER>(*this);
167 AffineMotionEstimationOptions & burtFilterStrength(
double strength)
169 vigra_precondition(0.25 <= strength && strength <= 0.5,
170 "AffineMotionEstimationOptions::burtFilterStrength(): strength must be between 0.25 and 0.5 (inclusive).");
171 burt_filter_strength = strength;
175 AffineMotionEstimationOptions & highestPyramidLevel(
unsigned int level)
177 highest_level = (int)level;
181 AffineMotionEstimationOptions & iterationsPerLevel(
unsigned int iter)
183 vigra_precondition(0 < iter,
184 "AffineMotionEstimationOptions::iterationsPerLevel(): must do at least one iteration per level.");
185 iterations_per_level = (int)iter;
189 AffineMotionEstimationOptions & useGaussianPyramid(
bool f =
true)
191 use_laplacian_pyramid = !f;
195 AffineMotionEstimationOptions & useLaplacianPyramid(
bool f =
true)
197 use_laplacian_pyramid = f;
204 struct TranslationEstimationFunctor
206 template <
class SplineImage,
class Image>
207 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const 209 int w = dest.width();
210 int h = dest.height();
212 Matrix<double> grad(2,1), m(2,2), r(2,1), s(2,1);
213 double dx = matrix(0,0), dy = matrix(1,0);
215 for(
int y = 0; y < h; ++y)
217 double sx = matrix(0,1)*y + matrix(0,2);
218 double sy = matrix(1,1)*y + matrix(1,2);
219 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
221 if(!src.isInside(sx, sy))
224 grad(0,0) = src.dx(sx, sy);
225 grad(1,0) = src.dy(sx, sy);
226 double diff = dest(x, y) - src(sx, sy);
235 matrix(0,2) -= s(0,0);
236 matrix(1,2) -= s(1,0);
240 struct SimilarityTransformEstimationFunctor
242 template <
class SplineImage,
class Image>
243 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const 245 int w = dest.width();
246 int h = dest.height();
248 Matrix<double> grad(2,1), coord(4, 2), c(4, 1), m(4, 4), r(4,1), s(4,1);
251 double dx = matrix(0,0), dy = matrix(1,0);
253 for(
int y = 0; y < h; ++y)
255 double sx = matrix(0,1)*y + matrix(0,2);
256 double sy = matrix(1,1)*y + matrix(1,2);
257 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
259 if(!src.isInside(sx, sy))
262 grad(0,0) = src.dx(sx, sy);
263 grad(1,0) = src.dy(sx, sy);
264 coord(2,0) = (double)x;
265 coord(3,1) = (double)x;
266 coord(3,0) = -(double)y;
267 coord(2,1) = (double)y;
268 double diff = dest(x, y) - src(sx, sy);
278 matrix(0,2) -= s(0,0);
279 matrix(1,2) -= s(1,0);
280 matrix(0,0) -= s(2,0);
281 matrix(1,1) -= s(2,0);
282 matrix(0,1) += s(3,0);
283 matrix(1,0) -= s(3,0);
287 struct AffineTransformEstimationFunctor
289 template <
class SplineImage,
class Image>
290 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const 292 int w = dest.width();
293 int h = dest.height();
295 Matrix<double> grad(2,1), coord(6, 2), c(6, 1), m(6,6), r(6,1), s(6,1);
298 double dx = matrix(0,0), dy = matrix(1,0);
300 for(
int y = 0; y < h; ++y)
302 double sx = matrix(0,1)*y + matrix(0,2);
303 double sy = matrix(1,1)*y + matrix(1,2);
304 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
306 if(!src.isInside(sx, sy))
309 grad(0,0) = src.dx(sx, sy);
310 grad(1,0) = src.dy(sx, sy);
311 coord(2,0) = (double)x;
312 coord(4,1) = (double)x;
313 coord(3,0) = (double)y;
314 coord(5,1) = (double)y;
315 double diff = dest(x, y) - src(sx, sy);
325 matrix(0,2) -= s(0,0);
326 matrix(1,2) -= s(1,0);
327 matrix(0,0) -= s(2,0);
328 matrix(0,1) -= s(3,0);
329 matrix(1,0) -= s(4,0);
330 matrix(1,1) -= s(5,0);
334 template <
class SrcIterator,
class SrcAccessor,
335 class DestIterator,
class DestAccessor,
336 int SPLINEORDER,
class Functor>
338 estimateAffineMotionImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
339 DestIterator dul, DestIterator dlr, DestAccessor dest,
340 Matrix<double> & affineMatrix,
341 AffineMotionEstimationOptions<SPLINEORDER>
const & options,
344 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote STmpType;
346 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote DTmpType;
349 int toplevel = options.highest_level;
353 if(options.use_laplacian_pyramid)
364 Matrix<double> currentMatrix(affineMatrix(2,2) == 0.0
365 ? identityMatrix<double>(3)
367 currentMatrix(0,2) /= std::pow(2.0, toplevel);
368 currentMatrix(1,2) /= std::pow(2.0, toplevel);
370 for(
int level = toplevel; level >= 0; --level)
374 for(
int iter = 0; iter < options.iterations_per_level; ++iter)
376 motionModel(sp, destPyramid[level], currentMatrix);
381 currentMatrix(0,2) *= 2.0;
382 currentMatrix(1,2) *= 2.0;
386 affineMatrix = currentMatrix;
453 template <
class SrcIterator,
class SrcAccessor,
454 class DestIterator,
class DestAccessor,
457 estimateTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
458 DestIterator dul, DestIterator dlr, DestAccessor dest,
459 Matrix<double> & affineMatrix,
460 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
462 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
463 options, detail::TranslationEstimationFunctor());
466 template <
class SrcIterator,
class SrcAccessor,
467 class DestIterator,
class DestAccessor>
469 estimateTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
470 DestIterator dul, DestIterator dlr, DestAccessor dest,
471 Matrix<double> & affineMatrix)
473 estimateTranslation(sul, slr, src, dul, dlr, dest,
474 affineMatrix, AffineMotionEstimationOptions<>());
477 template <
class SrcIterator,
class SrcAccessor,
478 class DestIterator,
class DestAccessor,
481 estimateTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
482 triple<DestIterator, DestIterator, DestAccessor> dest,
483 Matrix<double> & affineMatrix,
484 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
486 estimateTranslation(src.first, src.second, src.third, dest.first, dest.second, dest.third,
487 affineMatrix, options);
490 template <
class SrcIterator,
class SrcAccessor,
491 class DestIterator,
class DestAccessor>
493 estimateTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
494 triple<DestIterator, DestIterator, DestAccessor> dest,
495 Matrix<double> & affineMatrix)
497 estimateTranslation(src.first, src.second, src.third, dest.first, dest.second, dest.third,
498 affineMatrix, AffineMotionEstimationOptions<>());
501 template <
class T1,
class S1,
507 Matrix<double> & affineMatrix,
508 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
510 estimateTranslation(srcImageRange(src), destImageRange(dest),
511 affineMatrix, options);
514 template <
class T1,
class S1,
519 Matrix<double> & affineMatrix)
521 estimateTranslation(srcImageRange(src), destImageRange(dest),
522 affineMatrix, AffineMotionEstimationOptions<>());
589 template <
class SrcIterator,
class SrcAccessor,
590 class DestIterator,
class DestAccessor,
593 estimateSimilarityTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
594 DestIterator dul, DestIterator dlr, DestAccessor dest,
595 Matrix<double> & affineMatrix,
596 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
598 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
599 options, detail::SimilarityTransformEstimationFunctor());
602 template <
class SrcIterator,
class SrcAccessor,
603 class DestIterator,
class DestAccessor>
605 estimateSimilarityTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
606 DestIterator dul, DestIterator dlr, DestAccessor dest,
607 Matrix<double> & affineMatrix)
609 estimateSimilarityTransform(sul, slr, src, dul, dlr, dest,
610 affineMatrix, AffineMotionEstimationOptions<>());
613 template <
class SrcIterator,
class SrcAccessor,
614 class DestIterator,
class DestAccessor,
617 estimateSimilarityTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
618 triple<DestIterator, DestIterator, DestAccessor> dest,
619 Matrix<double> & affineMatrix,
620 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
622 estimateSimilarityTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
623 affineMatrix, options);
626 template <
class SrcIterator,
class SrcAccessor,
627 class DestIterator,
class DestAccessor>
629 estimateSimilarityTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
630 triple<DestIterator, DestIterator, DestAccessor> dest,
631 Matrix<double> & affineMatrix)
633 estimateSimilarityTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
634 affineMatrix, AffineMotionEstimationOptions<>());
637 template <
class T1,
class S1,
643 Matrix<double> & affineMatrix,
644 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
646 estimateSimilarityTransform(srcImageRange(src), destImageRange(dest),
647 affineMatrix, options);
650 template <
class T1,
class S1,
655 Matrix<double> & affineMatrix)
657 estimateSimilarityTransform(srcImageRange(src), destImageRange(dest),
658 affineMatrix, AffineMotionEstimationOptions<>());
725 template <
class SrcIterator,
class SrcAccessor,
726 class DestIterator,
class DestAccessor,
729 estimateAffineTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
730 DestIterator dul, DestIterator dlr, DestAccessor dest,
731 Matrix<double> & affineMatrix,
732 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
734 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
735 options, detail::AffineTransformEstimationFunctor());
738 template <
class SrcIterator,
class SrcAccessor,
739 class DestIterator,
class DestAccessor>
741 estimateAffineTransform(SrcIterator sul, SrcIterator slr, SrcAccessor src,
742 DestIterator dul, DestIterator dlr, DestAccessor dest,
743 Matrix<double> & affineMatrix)
745 estimateAffineTransform(sul, slr, src, dul, dlr, dest,
746 affineMatrix, AffineMotionEstimationOptions<>());
749 template <
class SrcIterator,
class SrcAccessor,
750 class DestIterator,
class DestAccessor,
753 estimateAffineTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
754 triple<DestIterator, DestIterator, DestAccessor> dest,
755 Matrix<double> & affineMatrix,
756 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
758 estimateAffineTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
759 affineMatrix, options);
762 template <
class SrcIterator,
class SrcAccessor,
763 class DestIterator,
class DestAccessor>
765 estimateAffineTransform(triple<SrcIterator, SrcIterator, SrcAccessor> src,
766 triple<DestIterator, DestIterator, DestAccessor> dest,
767 Matrix<double> & affineMatrix)
769 estimateAffineTransform(src.first, src.second, src.third, dest.first, dest.second, dest.third,
770 affineMatrix, AffineMotionEstimationOptions<>());
773 template <
class T1,
class S1,
779 Matrix<double> & affineMatrix,
780 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
782 vigra_precondition(src.
shape() == dest.
shape(),
783 "estimateAffineTransform(): shape mismatch between input and output.");
784 estimateAffineTransform(srcImageRange(src), destImageRange(dest),
785 affineMatrix, options);
788 template <
class T1,
class S1,
793 Matrix<double> & affineMatrix)
795 vigra_precondition(src.
shape() == dest.
shape(),
796 "estimateAffineTransform(): shape mismatch between input and output.");
797 estimateAffineTransform(srcImageRange(src), destImageRange(dest),
798 affineMatrix, AffineMotionEstimationOptions<>());
const difference_type & shape() const
Definition: multi_array.hxx:1551
Definition: accessor.hxx:43
void estimateSimilarityTransform(...)
Estimate the optical flow between two images according to a similarity transform model (e...
void estimateAffineTransform(...)
Estimate the optical flow between two images according to an affine transform model (e...
void estimateTranslation(...)
Estimate the optical flow between two images according to a translation model.
void pyramidReduceBurtFilter(...)
Two-fold down-sampling for image pyramid construction.
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1457
Fundamental class template for images.
Definition: basicimage.hxx:473
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:655
void pyramidReduceBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Create a Laplacian pyramid.
Definition: resampling_convolution.hxx:1154
Create a continuous view onto a discrete image using splines.
Definition: splineimageview.hxx:99
Class template for logarithmically tapering image pyramids.
Definition: imagecontainer.hxx:465
linalg::TemporaryMatrix< double > affineMatrix2DFromCorrespondingPoints(SrcIterator s, SrcIterator send, DestIterator d)
Create homogeneous matrix that maps corresponding points onto each other.
Definition: affine_registration.hxx:70