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" 46 #include "affinegeometry.hxx" 71 template <
class SrcIterator,
class DestIterator>
72 linalg::TemporaryMatrix<double>
77 linalg::TemporaryMatrix<double> ret(identityMatrix<double>(3));
81 ret(0,2) = (*d)[0] - (*s)[0];
82 ret(1,2) = (*d)[1] - (*s)[1];
86 Matrix<double> m(4,4), r(4,1), so(4,1);
88 for(
int k=0; k<size; ++k, ++s, ++d)
100 r(2*k+1,0) = (*d)[1];
104 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
115 Matrix<double> m(3,3), rx(3,1), sx(3,1), ry(3,1), sy(3,1), c(3,1);
117 for(
int k=0; k<size; ++k, ++s, ++d)
128 vigra_fail(
"affineMatrix2DFromCorrespondingPoints(): singular solution matrix.");
147 template <
int SPLINEORDER = 2>
151 double burt_filter_strength;
152 int highest_level, iterations_per_level;
153 bool use_laplacian_pyramid;
156 : burt_filter_strength(0.4),
158 iterations_per_level(4),
159 use_laplacian_pyramid(
false)
164 : burt_filter_strength(other.burt_filter_strength),
165 highest_level(other.highest_level),
166 iterations_per_level(other.iterations_per_level),
167 use_laplacian_pyramid(other.use_laplacian_pyramid)
180 template <
int NEWORDER>
199 vigra_precondition(0.25 <= center && center <= 0.5,
200 "AffineMotionEstimationOptions::burtFilterCenterStrength(): center must be between 0.25 and 0.5 (inclusive).");
201 burt_filter_strength = center;
214 highest_level = (int)level;
224 vigra_precondition(0 < iter,
225 "AffineMotionEstimationOptions::iterationsPerLevel(): must do at least one iteration per level.");
226 iterations_per_level = (int)iter;
239 use_laplacian_pyramid = !f;
252 use_laplacian_pyramid = f;
259 struct TranslationEstimationFunctor
261 template <
class SplineImage,
class Image>
262 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const 264 int w = dest.width();
265 int h = dest.height();
267 Matrix<double> grad(2,1), m(2,2), r(2,1), s(2,1);
268 double dx = matrix(0,0), dy = matrix(1,0);
270 for(
int y = 0; y < h; ++y)
272 double sx = matrix(0,1)*y + matrix(0,2);
273 double sy = matrix(1,1)*y + matrix(1,2);
274 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
276 if(!src.isInside(sx, sy))
279 grad(0,0) = src.dx(sx, sy);
280 grad(1,0) = src.dy(sx, sy);
281 double diff = dest(x, y) - src(sx, sy);
290 matrix(0,2) -= s(0,0);
291 matrix(1,2) -= s(1,0);
295 struct SimilarityTransformEstimationFunctor
297 template <
class SplineImage,
class Image>
298 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const 300 int w = dest.width();
301 int h = dest.height();
303 Matrix<double> grad(2,1), coord(4, 2), c(4, 1), m(4, 4), r(4,1), s(4,1);
306 double dx = matrix(0,0), dy = matrix(1,0);
308 for(
int y = 0; y < h; ++y)
310 double sx = matrix(0,1)*y + matrix(0,2);
311 double sy = matrix(1,1)*y + matrix(1,2);
312 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
314 if(!src.isInside(sx, sy))
317 grad(0,0) = src.dx(sx, sy);
318 grad(1,0) = src.dy(sx, sy);
319 coord(2,0) = (double)x;
320 coord(3,1) = (double)x;
321 coord(3,0) = -(double)y;
322 coord(2,1) = (double)y;
323 double diff = dest(x, y) - src(sx, sy);
333 matrix(0,2) -= s(0,0);
334 matrix(1,2) -= s(1,0);
335 matrix(0,0) -= s(2,0);
336 matrix(1,1) -= s(2,0);
337 matrix(0,1) += s(3,0);
338 matrix(1,0) -= s(3,0);
342 struct AffineTransformEstimationFunctor
344 template <
class SplineImage,
class Image>
345 void operator()(SplineImage
const & src, Image
const & dest, Matrix<double> & matrix)
const 347 int w = dest.width();
348 int h = dest.height();
350 Matrix<double> grad(2,1), coord(6, 2), c(6, 1), m(6,6), r(6,1), s(6,1);
353 double dx = matrix(0,0), dy = matrix(1,0);
355 for(
int y = 0; y < h; ++y)
357 double sx = matrix(0,1)*y + matrix(0,2);
358 double sy = matrix(1,1)*y + matrix(1,2);
359 for(
int x = 0; x < w; ++x, sx += dx, sy += dy)
361 if(!src.isInside(sx, sy))
364 grad(0,0) = src.dx(sx, sy);
365 grad(1,0) = src.dy(sx, sy);
366 coord(2,0) = (double)x;
367 coord(4,1) = (double)x;
368 coord(3,0) = (double)y;
369 coord(5,1) = (double)y;
370 double diff = dest(x, y) - src(sx, sy);
380 matrix(0,2) -= s(0,0);
381 matrix(1,2) -= s(1,0);
382 matrix(0,0) -= s(2,0);
383 matrix(0,1) -= s(3,0);
384 matrix(1,0) -= s(4,0);
385 matrix(1,1) -= s(5,0);
389 template <
class SrcIterator,
class SrcAccessor,
390 class DestIterator,
class DestAccessor,
391 int SPLINEORDER,
class Functor>
393 estimateAffineMotionImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
394 DestIterator dul, DestIterator dlr, DestAccessor dest,
395 Matrix<double> & affineMatrix,
399 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote STmpType;
401 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote DTmpType;
404 int toplevel = options.highest_level;
408 if(options.use_laplacian_pyramid)
419 Matrix<double> currentMatrix(affineMatrix(2,2) == 0.0
420 ? identityMatrix<double>(3)
422 currentMatrix(0,2) /= std::pow(2.0, toplevel);
423 currentMatrix(1,2) /= std::pow(2.0, toplevel);
425 for(
int level = toplevel; level >= 0; --level)
429 for(
int iter = 0; iter < options.iterations_per_level; ++iter)
431 motionModel(sp, destPyramid[level], currentMatrix);
436 currentMatrix(0,2) *= 2.0;
437 currentMatrix(1,2) *= 2.0;
441 affineMatrix = currentMatrix;
510 template <
class SrcIterator,
class SrcAccessor,
511 class DestIterator,
class DestAccessor,
515 DestIterator dul, DestIterator dlr, DestAccessor dest,
516 Matrix<double> & affineMatrix,
519 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
520 options, detail::TranslationEstimationFunctor());
523 template <
class SrcIterator,
class SrcAccessor,
524 class DestIterator,
class DestAccessor>
527 DestIterator dul, DestIterator dlr, DestAccessor dest,
528 Matrix<double> & affineMatrix)
534 template <
class SrcIterator,
class SrcAccessor,
535 class DestIterator,
class DestAccessor,
539 triple<DestIterator, DestIterator, DestAccessor> dest,
540 Matrix<double> & affineMatrix,
544 affineMatrix, options);
547 template <
class SrcIterator,
class SrcAccessor,
548 class DestIterator,
class DestAccessor>
551 triple<DestIterator, DestIterator, DestAccessor> dest,
552 Matrix<double> & affineMatrix)
558 template <
class T1,
class S1,
564 Matrix<double> & affineMatrix,
568 affineMatrix, options);
571 template <
class T1,
class S1,
576 Matrix<double> & affineMatrix)
648 template <
class SrcIterator,
class SrcAccessor,
649 class DestIterator,
class DestAccessor,
653 DestIterator dul, DestIterator dlr, DestAccessor dest,
654 Matrix<double> & affineMatrix,
657 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
658 options, detail::SimilarityTransformEstimationFunctor());
661 template <
class SrcIterator,
class SrcAccessor,
662 class DestIterator,
class DestAccessor>
665 DestIterator dul, DestIterator dlr, DestAccessor dest,
666 Matrix<double> & affineMatrix)
672 template <
class SrcIterator,
class SrcAccessor,
673 class DestIterator,
class DestAccessor,
677 triple<DestIterator, DestIterator, DestAccessor> dest,
678 Matrix<double> & affineMatrix,
682 affineMatrix, options);
685 template <
class SrcIterator,
class SrcAccessor,
686 class DestIterator,
class DestAccessor>
689 triple<DestIterator, DestIterator, DestAccessor> dest,
690 Matrix<double> & affineMatrix)
696 template <
class T1,
class S1,
702 Matrix<double> & affineMatrix,
706 affineMatrix, options);
709 template <
class T1,
class S1,
714 Matrix<double> & affineMatrix)
815 template <
class SrcIterator,
class SrcAccessor,
816 class DestIterator,
class DestAccessor,
820 DestIterator dul, DestIterator dlr, DestAccessor dest,
821 Matrix<double> & affineMatrix,
824 detail::estimateAffineMotionImpl(sul, slr, src, dul, dlr, dest, affineMatrix,
825 options, detail::AffineTransformEstimationFunctor());
828 template <
class SrcIterator,
class SrcAccessor,
829 class DestIterator,
class DestAccessor>
832 DestIterator dul, DestIterator dlr, DestAccessor dest,
833 Matrix<double> & affineMatrix)
839 template <
class SrcIterator,
class SrcAccessor,
840 class DestIterator,
class DestAccessor,
844 triple<DestIterator, DestIterator, DestAccessor> dest,
845 Matrix<double> & affineMatrix,
849 affineMatrix, options);
852 template <
class SrcIterator,
class SrcAccessor,
853 class DestIterator,
class DestAccessor>
856 triple<DestIterator, DestIterator, DestAccessor> dest,
857 Matrix<double> & affineMatrix)
863 template <
class T1,
class S1,
869 Matrix<double> & affineMatrix,
872 vigra_precondition(src.
shape() == dest.
shape(),
873 "estimateAffineTransform(): shape mismatch between input and output.");
875 affineMatrix, options);
878 template <
class T1,
class S1,
883 Matrix<double> & affineMatrix)
885 vigra_precondition(src.
shape() == dest.
shape(),
886 "estimateAffineTransform(): shape mismatch between input and output.");
AffineMotionEstimationOptions & highestPyramidLevel(unsigned int level)
Define the highest level of the image pyramid.
Definition: affine_registration.hxx:212
AffineMotionEstimationOptions & useGaussianPyramid(bool f=true)
Base registration on a Gaussian pyramid.
Definition: affine_registration.hxx:237
AffineMotionEstimationOptions< NEWORDER > splineOrder() const
Change the spline order for intensity interpolation.
Definition: affine_registration.hxx:181
Option object for affine registration functions.
Definition: affine_registration.hxx:148
AffineMotionEstimationOptions & iterationsPerLevel(unsigned int iter)
Number of Gauss-Newton iterations per level.
Definition: affine_registration.hxx:222
const difference_type & shape() const
Definition: multi_array.hxx:1596
AffineMotionEstimationOptions & burtFilterCenterStrength(double center)
Define amount of smoothing before subsampling to the next pyramid level.
Definition: affine_registration.hxx:197
Definition: accessor.hxx:43
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
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
AffineMotionEstimationOptions & useLaplacianPyramid(bool f=true)
Base registration on a Laplacian pyramid.
Definition: affine_registration.hxx:250
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
void pyramidReduceBurtLaplacian(ImagePyramid< Image, Alloc > &pyramid, int fromLevel, int toLevel, double centerValue=0.4)
Create a Laplacian pyramid.
Definition: resampling_convolution.hxx:1197
Create a continuous view onto a discrete image using splines.
Definition: splineimageview.hxx:99
Class template for logarithmically tapering image pyramids.
Definition: imagecontainer.hxx:468
linalg::TemporaryMatrix< double > affineMatrix2DFromCorrespondingPoints(SrcIterator s, SrcIterator send, DestIterator d)
Create homogeneous matrix that maps corresponding points onto each other.
Definition: affine_registration.hxx:73