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"
109 template <
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);
149 template <
class SplineImage,
150 class DestIterator,
class DestAccessor>
153 vigra::triple<DestIterator, DestIterator, DestAccessor> dest)
158 template <
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));
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);
333 template <
class SrcIterator,
class SrcAccessor,
334 class DestIterator,
class DestAccessor>
337 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
338 Matrix<double> & affineMatrix,
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;
401 template <
class SrcIterator,
class SrcAccessor,
402 class DestIterator,
class DestAccessor>
405 triple<DestIterator, DestIterator, DestAccessor> dest,
406 Matrix<double> & affineMatrix,
407 double & correlation_coefficient)
410 dest.first, dest.second, dest.third,
412 correlation_coefficient);
415 template <
class T1,
class S1,
420 Matrix<double> & affineMatrix,
421 double & correlation_coefficent)
424 destImageRange(dest),
426 correlation_coefficent);
502 template <
class SrcIterator,
class SrcAccessor,
503 class DestIterator,
class DestAccessor>
505 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
506 Matrix<double> & affine_matrix,
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];
581 template <
class SrcIterator,
class SrcAccessor,
582 class DestIterator,
class DestAccessor>
585 triple<DestIterator, DestIterator, DestAccessor> dest,
586 Matrix<double> & affineMatrix,
587 double & correlation_coefficent,
591 dest.first, dest.second, dest.third,
593 correlation_coefficent,
597 template <
class T1,
class S1,
602 Matrix<double> & affineMatrix,
603 double & correlation_coefficent,
607 destImageRange(dest),
609 correlation_coefficent,
677 template <
class SrcIterator,
class SrcAccessor,
678 class DestIterator,
class DestAccessor>
681 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
682 Matrix<double> & affineMatrix,
683 double & rotation_correlation,
684 double & translation_correlation,
687 typename DestIterator::difference_type d_shape = d_lr - d_ul;
690 Matrix<double> rotation_matrix;
692 srcIterRange(d_ul+border, d_lr-border, d_acc),
694 rotation_correlation);
702 Matrix<double> translation_matrix;
704 srcIterRange(d_ul, d_lr, d_acc),
706 translation_correlation,
709 affineMatrix = rotation_matrix * translation_matrix;
712 template <
class SrcIterator,
class SrcAccessor,
713 class DestIterator,
class DestAccessor>
716 triple<DestIterator, DestIterator, DestAccessor> dest,
717 Matrix<double> & affineMatrix,
718 double & rotation_correlation,
719 double & translation_correlation,
723 dest.first, dest.second, dest.third,
725 rotation_correlation,
726 translation_correlation,
730 template <
class T1,
class S1,
735 Matrix<double> & affineMatrix,
736 double & rotation_correlation,
737 double & translation_correlation,
741 destImageRange(dest),
743 rotation_correlation,
744 translation_correlation,
void fastNormalizedCrossCorrelation(...)
This function performes a fast normalized cross-correlation.
int y
Definition: diff2d.hxx:392
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
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
int x
Definition: diff2d.hxx:385
Two dimensional difference vector.
Definition: diff2d.hxx:185
void estimateGlobalRotationTranslation(...)
Estimate the (global) rotation and translation between two images by means a normalized cross correla...
void estimateGlobalRotation(...)
Estimate the rotation between two images by means of a normalized cross correlation matching of the F...
BasicImageView< T > makeBasicImageView(MultiArrayView< 2, T, Stride > const &array)
Definition: multi_array.hxx:3476
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void estimateGlobalTranslation(...)
Estimate the translation between two images by means of a normalized cross correlation matching...
void copyImage(...)
Copy source image into destination image.
Definition: fftw3.hxx:1377
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
Fundamental class template for images.
Definition: basicimage.hxx:475
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:704
Create a continuous view onto a discrete image using splines.
Definition: splineimageview.hxx:99
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
void normalizedCrossCorrelation(...)
This function performes a (slow) normalized cross-correlation.