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"
65 template <
class SrcIterator,
class DestIterator>
66 linalg::TemporaryMatrix<double>
71 linalg::TemporaryMatrix<double> ret(identityMatrix<double>(3));
75 ret(0,2) = (*d)[0] - (*s)[0];
76 ret(1,2) = (*d)[1] - (*s)[1];
80 Matrix<double> m(4,4), r(4,1), so(4,1);
82 for(
int k=0; k<size; ++k, ++s, ++d)
98 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
109 Matrix<double> m(3,3), rx(3,1), sx(3,1), ry(3,1), sy(3,1), c(3,1);
111 for(
int k=0; k<size; ++k, ++s, ++d)
122 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
135 template <
int SPLINEORDER = 2>
136 class AffineMotionEstimationOptions
139 double burt_filter_strength;
140 int highest_level, iterations_per_level;
141 bool use_laplacian_pyramid;
143 AffineMotionEstimationOptions()
144 : burt_filter_strength(0.4),
146 iterations_per_level(4),
147 use_laplacian_pyramid(false)
151 AffineMotionEstimationOptions(AffineMotionEstimationOptions<ORDER>
const & other)
152 : burt_filter_strength(other.burt_filter_strength),
153 highest_level(other.highest_level),
154 iterations_per_level(other.iterations_per_level),
155 use_laplacian_pyramid(other.use_laplacian_pyramid)
158 template <
int NEWORDER>
159 AffineMotionEstimationOptions<NEWORDER> splineOrder()
const
161 return AffineMotionEstimationOptions<NEWORDER>(*this);
164 AffineMotionEstimationOptions & burtFilterStrength(
double strength)
166 vigra_precondition(0.25 <= strength && strength <= 0.5,
167 "AffineMotionEstimationOptions::burtFilterStrength(): strength must be between 0.25 and 0.5 (inclusive).");
168 burt_filter_strength = strength;
172 AffineMotionEstimationOptions & highestPyramidLevel(
unsigned int level)
174 highest_level = (int)level;
178 AffineMotionEstimationOptions & iterationsPerLevel(
unsigned int iter)
180 vigra_precondition(0 < iter,
181 "AffineMotionEstimationOptions::iterationsPerLevel(): must do at least one iteration per level.");
182 iterations_per_level = (int)iter;
186 AffineMotionEstimationOptions & useGaussianPyramid(
bool f =
true)
188 use_laplacian_pyramid = !f;
192 AffineMotionEstimationOptions & useLaplacianPyramid(
bool f =
true)
194 use_laplacian_pyramid = f;
201 struct TranslationEstimationFunctor
203 template <
class SplineImage,
class Image>
204 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const
206 int w = dest.width();
207 int h = dest.height();
209 Matrix<double> grad(2,1), m(2,2), r(2,1), s(2,1);
210 double dx = matrix(0,0), dy = matrix(1,0);
212 for(
int y = 0; y < h; ++y)
214 double sx = matrix(0,1)*y + matrix(0,2);
215 double sy = matrix(1,1)*y + matrix(1,2);
216 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
218 if(!src.isInside(sx, sy))
221 grad(0,0) = src.dx(sx, sy);
222 grad(1,0) = src.dy(sx, sy);
223 double diff = dest(x, y) - src(sx, sy);
232 matrix(0,2) -= s(0,0);
233 matrix(1,2) -= s(1,0);
237 struct SimilarityTransformEstimationFunctor
239 template <
class SplineImage,
class Image>
240 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const
242 int w = dest.width();
243 int h = dest.height();
245 Matrix<double> grad(2,1), coord(4, 2), c(4, 1), m(4, 4), r(4,1), s(4,1);
248 double dx = matrix(0,0), dy = matrix(1,0);
250 for(
int y = 0; y < h; ++y)
252 double sx = matrix(0,1)*y + matrix(0,2);
253 double sy = matrix(1,1)*y + matrix(1,2);
254 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
256 if(!src.isInside(sx, sy))
259 grad(0,0) = src.dx(sx, sy);
260 grad(1,0) = src.dy(sx, sy);
261 coord(2,0) = (double)x;
262 coord(3,1) = (double)x;
263 coord(3,0) = -(double)y;
264 coord(2,1) = (double)y;
265 double diff = dest(x, y) - src(sx, sy);
275 matrix(0,2) -= s(0,0);
276 matrix(1,2) -= s(1,0);
277 matrix(0,0) -= s(2,0);
278 matrix(1,1) -= s(2,0);
279 matrix(0,1) += s(3,0);
280 matrix(1,0) -= s(3,0);
284 struct AffineTransformEstimationFunctor
286 template <
class SplineImage,
class Image>
287 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const
289 int w = dest.width();
290 int h = dest.height();
292 Matrix<double> grad(2,1), coord(6, 2), c(6, 1), m(6,6), r(6,1), s(6,1);
295 double dx = matrix(0,0), dy = matrix(1,0);
297 for(
int y = 0; y < h; ++y)
299 double sx = matrix(0,1)*y + matrix(0,2);
300 double sy = matrix(1,1)*y + matrix(1,2);
301 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
303 if(!src.isInside(sx, sy))
306 grad(0,0) = src.dx(sx, sy);
307 grad(1,0) = src.dy(sx, sy);
308 coord(2,0) = (double)x;
309 coord(4,1) = (double)x;
310 coord(3,0) = (double)y;
311 coord(5,1) = (double)y;
312 double diff = dest(x, y) - src(sx, sy);
322 matrix(0,2) -= s(0,0);
323 matrix(1,2) -= s(1,0);
324 matrix(0,0) -= s(2,0);
325 matrix(0,1) -= s(3,0);
326 matrix(1,0) -= s(4,0);
327 matrix(1,1) -= s(5,0);
331 template <
class SrcIterator,
class SrcAccessor,
332 class DestIterator,
class DestAccessor,
333 int SPLINEORDER,
class Functor>
335 estimateAffineMotionImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
336 DestIterator dul, DestIterator dlr, DestAccessor dest,
337 Matrix<double> & affineMatrix,
338 AffineMotionEstimationOptions<SPLINEORDER>
const & options,
341 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote STmpType;
343 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote DTmpType;
346 int toplevel = options.highest_level;
350 if(options.use_laplacian_pyramid)
361 Matrix<double> currentMatrix(affineMatrix(2,2) == 0.0
362 ? identityMatrix<double>(3)
364 currentMatrix(0,2) /=
std::pow(2.0, toplevel);
365 currentMatrix(1,2) /=
std::pow(2.0, toplevel);
367 for(
int level = toplevel; level >= 0; --level)
371 for(
int iter = 0; iter < options.iterations_per_level; ++iter)
373 motionModel(sp, destPyramid[level], currentMatrix);
378 currentMatrix(0,2) *= 2.0;
379 currentMatrix(1,2) *= 2.0;
383 affineMatrix = currentMatrix;
436 template <
class SrcIterator,
class SrcAccessor,
437 class DestIterator,
class DestAccessor,
441 DestIterator dul, DestIterator dlr, DestAccessor dest,
442 Matrix<double> & affineMatrix,
443 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
445 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
446 options, detail::TranslationEstimationFunctor());
449 template <
class SrcIterator,
class SrcAccessor,
450 class DestIterator,
class DestAccessor>
453 DestIterator dul, DestIterator dlr, DestAccessor dest,
454 Matrix<double> & affineMatrix)
457 affineMatrix, AffineMotionEstimationOptions<>());
460 template <
class SrcIterator,
class SrcAccessor,
461 class DestIterator,
class DestAccessor,
465 triple<DestIterator, DestIterator, DestAccessor> dest,
466 Matrix<double> & affineMatrix,
467 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
470 affineMatrix, options);
473 template <
class SrcIterator,
class SrcAccessor,
474 class DestIterator,
class DestAccessor>
477 triple<DestIterator, DestIterator, DestAccessor> dest,
478 Matrix<double> & affineMatrix)
481 affineMatrix, AffineMotionEstimationOptions<>());
533 template <
class SrcIterator,
class SrcAccessor,
534 class DestIterator,
class DestAccessor,
538 DestIterator dul, DestIterator dlr, DestAccessor dest,
539 Matrix<double> & affineMatrix,
540 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
542 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
543 options, detail::SimilarityTransformEstimationFunctor());
546 template <
class SrcIterator,
class SrcAccessor,
547 class DestIterator,
class DestAccessor>
550 DestIterator dul, DestIterator dlr, DestAccessor dest,
551 Matrix<double> & affineMatrix)
554 affineMatrix, AffineMotionEstimationOptions<>());
557 template <
class SrcIterator,
class SrcAccessor,
558 class DestIterator,
class DestAccessor,
562 triple<DestIterator, DestIterator, DestAccessor> dest,
563 Matrix<double> & affineMatrix,
564 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
567 affineMatrix, options);
570 template <
class SrcIterator,
class SrcAccessor,
571 class DestIterator,
class DestAccessor>
574 triple<DestIterator, DestIterator, DestAccessor> dest,
575 Matrix<double> & affineMatrix)
578 affineMatrix, AffineMotionEstimationOptions<>());
630 template <
class SrcIterator,
class SrcAccessor,
631 class DestIterator,
class DestAccessor,
635 DestIterator dul, DestIterator dlr, DestAccessor dest,
636 Matrix<double> & affineMatrix,
637 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
639 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
640 options, detail::AffineTransformEstimationFunctor());
643 template <
class SrcIterator,
class SrcAccessor,
644 class DestIterator,
class DestAccessor>
647 DestIterator dul, DestIterator dlr, DestAccessor dest,
648 Matrix<double> & affineMatrix)
651 affineMatrix, AffineMotionEstimationOptions<>());
654 template <
class SrcIterator,
class SrcAccessor,
655 class DestIterator,
class DestAccessor,
659 triple<DestIterator, DestIterator, DestAccessor> dest,
660 Matrix<double> & affineMatrix,
661 AffineMotionEstimationOptions<SPLINEORDER>
const & options)
664 affineMatrix, options);
667 template <
class SrcIterator,
class SrcAccessor,
668 class DestIterator,
class DestAccessor>
671 triple<DestIterator, DestIterator, DestAccessor> dest,
672 Matrix<double> & affineMatrix)
675 affineMatrix, AffineMotionEstimationOptions<>());