36#ifndef VIGRA_AFFINE_REGISTRATION_FFT_HXX
37#define VIGRA_AFFINE_REGISTRATION_FFT_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"
47#include "correlation.hxx"
109template <
class SplineImage,
110 class DestIterator,
class DestAccessor>
113 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc)
115 typename DestIterator::difference_type d_shape = (d_lr - d_ul);
117 int s_w = src.width(),
120 int s_size = min(s_w, s_h);
125 double r_max = s_size / 2.0;
127 DestIterator yd = d_ul;
128 DestIterator xd = yd;
130 for (
int t_step = 0; t_step < d_h; t_step++, yd.y++)
133 for (
int r_step = 0; r_step < d_w; r_step++, xd.x++)
135 double theta = 2.0 * M_PI * double(t_step) / double(d_h);
136 double r = r_max * double(r_step) / double(d_w);
137 double u = r * cos(theta) + r_max;
138 double v = r * -sin(theta) + r_max;
140 if ( u >= 0 && u < s_size
141 && v >= 0 && v < s_size)
143 d_acc.set(src(u, v), xd);
149template <
class SplineImage,
150 class DestIterator,
class DestAccessor>
153 vigra::triple<DestIterator, DestIterator, DestAccessor> dest)
158template <
class SplineImage,
162 MultiArrayView<2, T1, S1> dest)
172 template <
class SrcIterator,
class SrcAccessor,
173 class DestIterator,
class DestAccessor>
175 maximumFastNCC(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
176 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
177 TinyVector<int,2> & position,
178 double & correlation_coefficent)
180 typename DestIterator::difference_type s_shape = s_lr - s_ul;
181 typename DestIterator::difference_type d_shape = d_lr - d_ul;
183 MultiArray<2, float> src(s_shape.x, s_shape.y), dest(d_shape.x, d_shape.y), ncc(d_shape.x, d_shape.y);
187 copyImage(srcIterRange(s_ul, s_lr, s_acc), destImage(src_view));
188 copyImage(srcIterRange(d_ul, d_lr, d_acc), destImage(dest_view));
195 float max_val = -1.0;
197 for (
int y = 0; y < ncc.height()-s_shape.y; y++)
199 for (
int x = 0; x < ncc.width()-s_shape.x; x++)
201 val = ncc(x+s_shape.x/2, y+s_shape.y/2);
215 correlation_coefficent = max_val;
218 template <
class SrcIterator,
class SrcAccessor,
219 class DestIterator,
class DestAccessor>
221 maximumFastNCC(triple<SrcIterator, SrcIterator, SrcAccessor> src,
222 triple<DestIterator, DestIterator, DestAccessor> dest,
223 TinyVector<int,2> & position,
224 double & correlation_coefficent)
226 maximumFastNCC(src.first, src.second, src.third,
227 dest.first, dest.second, dest.third,
229 correlation_coefficent);
232 template <
class SrcIterator,
class SrcAccessor,
233 class DestIterator,
class DestAccessor>
234 void fourierLogAbsSpectrumInPolarCoordinates(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
235 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc)
237 using namespace vigra;
239 typename SrcIterator::difference_type shape = s_lr - s_ul;
250 destIterRange(d_ul, d_lr, d_acc));
252 transformImage(srcIterRange(d_ul,d_lr,d_acc), destIter(d_ul,d_acc), log(
abs(functor::Arg1())));
255 template <
class SrcIterator,
class SrcAccessor,
256 class DestIterator,
class DestAccessor>
257 void fourierLogAbsSpectrumInPolarCoordinates(triple<SrcIterator, SrcIterator, SrcAccessor> src,
258 triple<DestIterator, DestIterator, DestAccessor> dest)
260 fourierLogAbsSpectrumInPolarCoordinates(src.first, src.second, src.third, dest.first, dest.second, dest.third);
333template <
class SrcIterator,
class SrcAccessor,
334 class DestIterator,
class DestAccessor>
337 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
339 double & correlation_coefficient)
342 typename SrcIterator::difference_type s_shape = s_lr - s_ul;
343 Diff2D s2_shape(min(s_shape.x, s_shape.y),min(s_shape.x, s_shape.y));
344 Diff2D s2_offset = (s_shape-s2_shape)/2;
346 typename DestIterator::difference_type d_shape = d_lr - d_ul;
347 Diff2D d2_shape(min(d_shape.x, d_shape.y),min(d_shape.x, d_shape.y));
348 Diff2D d2_offset = (d_shape-d2_shape)/2;
351 Diff2D mean_shape = (s_shape + d_shape)/2;
353 int size = min(mean_shape.
x, mean_shape.
y);
357 const int pc_w = size,
361 DImage s_polCoords(pc_w, pc_h/2),
362 d_polCoords(pc_w, pc_h),
365 detail::fourierLogAbsSpectrumInPolarCoordinates(srcIterRange(s_ul+s2_offset, s_ul+s2_offset+s2_shape, s_acc),
366 destImageRange(d_polCoords));
367 copyImage(srcIterRange(d_polCoords.upperLeft(), d_polCoords.upperLeft() +
vigra::Diff2D(pc_w, pc_h/2), d_polCoords.accessor()),
368 destImage(s_polCoords));
370 detail::fourierLogAbsSpectrumInPolarCoordinates(srcIterRange(d_ul+d2_offset, d_ul+d2_offset+d2_shape, d_acc),
371 destImageRange(d_polCoords));
384 for (
int y=0; y<pc_h/2; y++)
386 val = ncc(x,y+pc_h/4);
395 double theta = double(max_idx) / pc_h * 360.0;
398 correlation_coefficient = max_val;
401template <
class SrcIterator,
class SrcAccessor,
402 class DestIterator,
class DestAccessor>
405 triple<DestIterator, DestIterator, DestAccessor> dest,
407 double & correlation_coefficient)
410 dest.first, dest.second, dest.third,
412 correlation_coefficient);
415template <
class T1,
class S1,
421 double & correlation_coefficent)
424 destImageRange(dest),
426 correlation_coefficent);
502template <
class SrcIterator,
class SrcAccessor,
503 class DestIterator,
class DestAccessor>
505 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
507 double & correlation_coefficent,
510 ignore_argument(d_lr);
511 typename SrcIterator::difference_type s_shape = s_lr - s_ul;
514 Diff2D q_shape = (s_shape - border - border)/2;
515 if (q_shape.
x % 2 == 0) q_shape.
x--;
516 if (q_shape.
y % 2 == 0) q_shape.
y--;
519 q_offsets[0] = border;
520 q_offsets[1] =
Diff2D(s_shape.x, 0)/2 + border;
521 q_offsets[2] = s_shape/4;
522 q_offsets[3] =
Diff2D(0, s_shape.y)/2 + border;
523 q_offsets[4] = s_shape/2 + border;
526 double translation_correlations[5] = {0.0,0.0,0.0,0.0,0.0};
527 int translation_votes[5] = {1,1,1,1,1};
531 for (
int q=0; q!=5; q++)
533 Diff2D offset = q_offsets[q];
536 detail::maximumFastNCC(srcIterRange(d_ul+offset, d_ul+offset+q_shape, d_acc),
537 srcIterRange(s_ul, s_lr, s_acc),
538 translation_vectors[q],
539 translation_correlations[q]);
543 if(translation_correlations[q] > translation_correlations[max_corr_idx])
548 for (
int q_old=0; q_old!=q; q_old++)
550 if (translation_vectors[q] == translation_vectors[q_old])
552 translation_votes[q_old]++;
559 for (
int q=0; q!=5; q++)
561 if(translation_votes[q] > translation_votes[max_votes_idx])
568 if(translation_votes[max_votes_idx] > 1)
570 best_idx = max_votes_idx;
574 best_idx = max_corr_idx;
578 correlation_coefficent = translation_correlations[best_idx];
581template <
class SrcIterator,
class SrcAccessor,
582 class DestIterator,
class DestAccessor>
585 triple<DestIterator, DestIterator, DestAccessor> dest,
587 double & correlation_coefficent,
591 dest.first, dest.second, dest.third,
593 correlation_coefficent,
597template <
class T1,
class S1,
603 double & correlation_coefficent,
607 destImageRange(dest),
609 correlation_coefficent,
677template <
class SrcIterator,
class SrcAccessor,
678 class DestIterator,
class DestAccessor>
681 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
683 double & rotation_correlation,
684 double & translation_correlation,
687 typename DestIterator::difference_type d_shape = d_lr - d_ul;
692 srcIterRange(d_ul+border, d_lr-border, d_acc),
694 rotation_correlation);
704 srcIterRange(d_ul, d_lr, d_acc),
706 translation_correlation,
709 affineMatrix = rotation_matrix * translation_matrix;
712template <
class SrcIterator,
class SrcAccessor,
713 class DestIterator,
class DestAccessor>
716 triple<DestIterator, DestIterator, DestAccessor> dest,
718 double & rotation_correlation,
719 double & translation_correlation,
723 dest.first, dest.second, dest.third,
725 rotation_correlation,
726 translation_correlation,
730template <
class T1,
class S1,
736 double & rotation_correlation,
737 double & translation_correlation,
741 destImageRange(dest),
743 rotation_correlation,
744 translation_correlation,
Fundamental class template for images.
Definition basicimage.hxx:476
Two dimensional difference vector.
Definition diff2d.hxx:186
int y
Definition diff2d.hxx:392
int x
Definition diff2d.hxx:385
Definition fftw3.hxx:1378
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
Create a continuous view onto a discrete image using splines.
Definition splineimageview.hxx:100
Class for fixed size vectors.
Definition tinyvector.hxx:1008
Definition matrix.hxx:125
void fourierTransform(...)
Compute forward and inverse Fourier transforms.
BasicImageView< T > makeBasicImageView(MultiArrayView< 2, T, Stride > const &array)
Definition multi_array.hxx:3478
void fastNormalizedCrossCorrelation(...)
This function performes a fast normalized cross-correlation.
void estimateGlobalTranslation(...)
Estimate the translation between two images by means of a normalized cross correlation matching.
linalg::TemporaryMatrix< double > rotationMatrix2DDegrees(double angle)
Create homogeneous matrix representing a 2D rotation about the coordinate origin.
Definition affinegeometry.hxx:133
linalg::TemporaryMatrix< double > translationMatrix2D(TinyVector< double, 2 > const &shift)
Create homogeneous matrix representing a 2D translation.
Definition affinegeometry.hxx:64
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
void moveDCToCenter(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image center.
void affineWarpImage(...)
Warp an image according to an affine transformation.
void transformToPolarCoordinates(SplineImage const &src, DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc)
Transforms a given image to its (image-centered) polar coordinates representation.
Definition affine_registration_fft.hxx:112
void estimateGlobalRotationTranslation(...)
Estimate the (global) rotation and translation between two images by means a normalized cross correla...
void normalizedCrossCorrelation(...)
This function performes a (slow) normalized cross-correlation.
void estimateGlobalRotation(...)
Estimate the rotation between two images by means of a normalized cross correlation matching of the F...
void copyImage(...)
Copy source image into destination image.