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;
345 typedef BasicImage<STmpType> STmpImage;
346 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote DTmpType;
347 typedef BasicImage<DTmpType> DTmpImage;
349 int toplevel = options.highest_level;
350 ImagePyramid<STmpImage> srcPyramid(0, toplevel, sul, slr, src);
351 ImagePyramid<DTmpImage> destPyramid(0, toplevel, dul, dlr, dest);
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)
372 SplineImageView<SPLINEORDER, STmpType> sp(srcImageRange(srcPyramid[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,
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>
470 DestIterator dul, DestIterator dlr, DestAccessor dest,
471 Matrix<double> & affineMatrix)
474 affineMatrix, AffineMotionEstimationOptions<>());
477 template <
class SrcIterator,
class SrcAccessor,
478 class DestIterator,
class DestAccessor,
482 triple<DestIterator, DestIterator, DestAccessor> dest,
483 Matrix<double> & affineMatrix,
484 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
487 affineMatrix, options);
490 template <
class SrcIterator,
class SrcAccessor,
491 class DestIterator,
class DestAccessor>
494 triple<DestIterator, DestIterator, DestAccessor> dest,
495 Matrix<double> & affineMatrix)
498 affineMatrix, AffineMotionEstimationOptions<>());
501 template <
class T1,
class S1,
506 MultiArrayView<2, T2, S2> dest,
507 Matrix<double> & affineMatrix,
508 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
511 affineMatrix, options);
514 template <
class T1,
class S1,
518 MultiArrayView<2, T2, S2> dest,
519 Matrix<double> & affineMatrix)
522 affineMatrix, AffineMotionEstimationOptions<>());
589 template <
class SrcIterator,
class SrcAccessor,
590 class DestIterator,
class DestAccessor,
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>
606 DestIterator dul, DestIterator dlr, DestAccessor dest,
607 Matrix<double> & affineMatrix)
610 affineMatrix, AffineMotionEstimationOptions<>());
613 template <
class SrcIterator,
class SrcAccessor,
614 class DestIterator,
class DestAccessor,
618 triple<DestIterator, DestIterator, DestAccessor> dest,
619 Matrix<double> & affineMatrix,
620 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
623 affineMatrix, options);
626 template <
class SrcIterator,
class SrcAccessor,
627 class DestIterator,
class DestAccessor>
630 triple<DestIterator, DestIterator, DestAccessor> dest,
631 Matrix<double> & affineMatrix)
634 affineMatrix, AffineMotionEstimationOptions<>());
637 template <
class T1,
class S1,
642 MultiArrayView<2, T2, S2> dest,
643 Matrix<double> & affineMatrix,
644 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
647 affineMatrix, options);
650 template <
class T1,
class S1,
654 MultiArrayView<2, T2, S2> dest,
655 Matrix<double> & affineMatrix)
658 affineMatrix, AffineMotionEstimationOptions<>());
725 template <
class SrcIterator,
class SrcAccessor,
726 class DestIterator,
class DestAccessor,
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>
742 DestIterator dul, DestIterator dlr, DestAccessor dest,
743 Matrix<double> & affineMatrix)
746 affineMatrix, AffineMotionEstimationOptions<>());
749 template <
class SrcIterator,
class SrcAccessor,
750 class DestIterator,
class DestAccessor,
754 triple<DestIterator, DestIterator, DestAccessor> dest,
755 Matrix<double> & affineMatrix,
756 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
759 affineMatrix, options);
762 template <
class SrcIterator,
class SrcAccessor,
763 class DestIterator,
class DestAccessor>
766 triple<DestIterator, DestIterator, DestAccessor> dest,
767 Matrix<double> & affineMatrix)
770 affineMatrix, AffineMotionEstimationOptions<>());
773 template <
class T1,
class S1,
778 MultiArrayView<2, T2, S2> dest,
779 Matrix<double> & affineMatrix,
780 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
782 vigra_precondition(src.shape() == dest.shape(),
783 "estimateAffineTransform(): shape mismatch between input and output.");
785 affineMatrix, options);
788 template <
class T1,
class S1,
792 MultiArrayView<2, T2, S2> dest,
793 Matrix<double> & affineMatrix)
795 vigra_precondition(src.shape() == dest.shape(),
796 "estimateAffineTransform(): shape mismatch between input and output.");
798 affineMatrix, AffineMotionEstimationOptions<>());