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;
342 typedef BasicImage<STmpType> STmpImage;
343 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote DTmpType;
344 typedef BasicImage<DTmpType> DTmpImage;
346 int toplevel = options.highest_level;
347 ImagePyramid<STmpImage> srcPyramid(0, toplevel, sul, slr, src);
348 ImagePyramid<DTmpImage> destPyramid(0, toplevel, dul, dlr, dest);
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)
369 SplineImageView<SPLINEORDER, STmpType> sp(srcImageRange(srcPyramid[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<>());
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
void pyramidReduceBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Create a Laplacian pyramid.
Definition: resampling_convolution.hxx:1012
linalg::TemporaryMatrix< double > affineMatrix2DFromCorrespondingPoints(SrcIterator s, SrcIterator send, DestIterator d)
Create homogeneous matrix that maps corresponding points onto each other.
Definition: affine_registration.hxx:67
bool linearSolve(const MultiArrayView< 2, T, C1 > &A, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &res, std::string method="QR")
Definition: linear_solve.hxx:1173