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;
451 doxygen_overloaded_function(template <...>
void estimateTranslation)
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<>());
587 doxygen_overloaded_function(template <...>
void estimateSimilarityTransform)
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<>());
723 doxygen_overloaded_function(template <...>
void estimateAffineTransform)
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<>());
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
void pyramidReduceBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Create a Laplacian pyramid.
Definition: resampling_convolution.hxx:1154
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