1 /************************************************************************/
2 /* */
3 /* Copyright 2007-2014 by Benjamin Seppke */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35
36 #ifndef VIGRA_AFFINE_REGISTRATION_FFT_HXX
37 #define VIGRA_AFFINE_REGISTRATION_FFT_HXX
38
39 #include "mathutil.hxx"
40 #include "matrix.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"
48
49 #include <cmath>
50
51 namespace vigra {
52
53 /** \addtogroup Registration Image Registration
54 */
55 //@{
56
57 /********************************************************/
58 /* */
59 /* transformToPolarCoordinates */
60 /* */
61 /********************************************************/
62
63 /** \brief Transforms a given image to its (image-centered) polar coordinates representation.
64
65 This algorithm transforms a given image (by means of an spline image view) to its
66 image-centered polar coordinates reprensentation. The sampling of the polar coordinate system
67 is determined by the shape of the dest. image.
68
69 <b> Declarations:</b>
70
71 <b>\#include</b> \<vigra/affine_registration_fft.hxx\><br>
72 Namespace: vigra
73
74 pass 2D array views:
75 \code
76 namespace vigra {
77 template <class SplineImage,
78 class T1, class S1>
79 void
80 transformToPolarCoordinates(SplineImage const & src,
81 MultiArrayView<2, T1, S1> dest);
82 }
83 \endcode
84
85 \deprecatedAPI{estimateTranslation}
86 pass \ref ImageIterators and \ref DataAccessors :
87 \code
88 namespace vigra {
89 template <class SplineImage,
90 class DestIterator, class DestAccessor>
91 void
92 transformToPolarCoordinates(SplineImage const & src,
93 DestIterator dul, DestIterator dlr, DestAccessor dest)
94 }
95 \endcode
96 use argument objects in conjunction with \ref ArgumentObjectFactories :
97 \code
98 namespace vigra {
99 template <class SplineImage,
100 class DestIterator, class DestAccessor>
101 void
102 transformToPolarCoordinates(SplineImage const & src,
103 triple<DestIterator, DestIterator, DestAccessor> dest)
104 }
105 \endcode
106 \deprecatedEnd
107 */
108
109 template <class SplineImage,
110 class DestIterator, class DestAccessor>
111 void
transformToPolarCoordinates(SplineImage const & src,DestIterator d_ul,DestIterator d_lr,DestAccessor d_acc)112 transformToPolarCoordinates(SplineImage const & src,
113 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc)
114 {
115 typename DestIterator::difference_type d_shape = (d_lr - d_ul);
116
117 int s_w = src.width(),
118 s_h = src.height();
119
120 int s_size = min(s_w, s_h);
121
122 int d_w = d_shape.x,
123 d_h = d_shape.y;
124
125 double r_max = s_size / 2.0;
126
127 DestIterator yd = d_ul;
128 DestIterator xd = yd;
129
130 for (int t_step = 0; t_step < d_h; t_step++, yd.y++)
131 {
132 xd = yd;
133 for (int r_step = 0; r_step < d_w; r_step++, xd.x++)
134 {
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;
139
140 if ( u >= 0 && u < s_size
141 && v >= 0 && v < s_size)
142 {
143 d_acc.set(src(u, v), xd);
144 }
145 }
146 }
147 }
148
149 template <class SplineImage,
150 class DestIterator, class DestAccessor>
151 inline void
transformToPolarCoordinates(SplineImage const & src,vigra::triple<DestIterator,DestIterator,DestAccessor> dest)152 transformToPolarCoordinates(SplineImage const & src,
153 vigra::triple<DestIterator, DestIterator, DestAccessor> dest)
154 {
155 transformToPolarCoordinates(src, dest.first, dest.second, dest.third);
156 }
157
158 template <class SplineImage,
159 class T1, class S1>
160 void
transformToPolarCoordinates(SplineImage const & src,MultiArrayView<2,T1,S1> dest)161 transformToPolarCoordinates(SplineImage const & src,
162 MultiArrayView<2, T1, S1> dest)
163 {
164 transformToPolarCoordinates(src, srcImageRange(dest));
165 }
166
167
168
169
170 namespace detail
171 {
172 template <class SrcIterator, class SrcAccessor,
173 class DestIterator, class DestAccessor>
174 void
maximumFastNCC(SrcIterator s_ul,SrcIterator s_lr,SrcAccessor s_acc,DestIterator d_ul,DestIterator d_lr,DestAccessor d_acc,TinyVector<int,2> & position,double & correlation_coefficent)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)
179 {
180 typename DestIterator::difference_type s_shape = s_lr - s_ul;
181 typename DestIterator::difference_type d_shape = d_lr - d_ul;
182
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);
184 BasicImageView<float> src_view = makeBasicImageView(src);
185 BasicImageView<float> dest_view = makeBasicImageView(dest);
186
187 copyImage(srcIterRange(s_ul, s_lr, s_acc), destImage(src_view));
188 copyImage(srcIterRange(d_ul, d_lr, d_acc), destImage(dest_view));
189
190 fastNormalizedCrossCorrelation(dest, src, ncc);
191
192 int max_x = 0;
193 int max_y = 0;
194 float val = 0.0;
195 float max_val = -1.0;
196
197 for (int y = 0; y < ncc.height()-s_shape.y; y++)
198 {
199 for (int x = 0; x < ncc.width()-s_shape.x; x++)
200 {
201 val = ncc(x+s_shape.x/2, y+s_shape.y/2);
202
203 if (val > max_val)
204 {
205 max_x = x;
206 max_y = y;
207 max_val = val;
208 }
209 }
210 }
211
212 position[0] = max_x;
213 position[1] = max_y;
214
215 correlation_coefficent = max_val;
216 }
217
218 template <class SrcIterator, class SrcAccessor,
219 class DestIterator, class DestAccessor>
220 inline void
maximumFastNCC(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,TinyVector<int,2> & position,double & correlation_coefficent)221 maximumFastNCC(triple<SrcIterator, SrcIterator, SrcAccessor> src,
222 triple<DestIterator, DestIterator, DestAccessor> dest,
223 TinyVector<int,2> & position,
224 double & correlation_coefficent)
225 {
226 maximumFastNCC(src.first, src.second, src.third,
227 dest.first, dest.second, dest.third,
228 position,
229 correlation_coefficent);
230 }
231
232 template <class SrcIterator, class SrcAccessor,
233 class DestIterator, class DestAccessor>
fourierLogAbsSpectrumInPolarCoordinates(SrcIterator s_ul,SrcIterator s_lr,SrcAccessor s_acc,DestIterator d_ul,DestIterator d_lr,DestAccessor d_acc)234 void fourierLogAbsSpectrumInPolarCoordinates(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
235 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc)
236 {
237 using namespace vigra;
238
239 typename SrcIterator::difference_type shape = s_lr - s_ul;
240
241 FFTWComplexImage fourier(shape);
242 FImage fft_mag(shape);
243
244 fourierTransform(srcIterRange(s_ul, s_lr, s_acc), destImage(fourier));
245 moveDCToCenter(srcImageRange(fourier, FFTWMagnitudeAccessor<>()), destImage(fft_mag));
246
247 vigra::SplineImageView<2, double> spl(srcImageRange(fft_mag));
248
249 transformToPolarCoordinates(spl,
250 destIterRange(d_ul, d_lr, d_acc));
251
252 transformImage(srcIterRange(d_ul,d_lr,d_acc), destIter(d_ul,d_acc), log(abs(functor::Arg1())));
253 }
254
255 template <class SrcIterator, class SrcAccessor,
256 class DestIterator, class DestAccessor>
fourierLogAbsSpectrumInPolarCoordinates(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest)257 void fourierLogAbsSpectrumInPolarCoordinates(triple<SrcIterator, SrcIterator, SrcAccessor> src,
258 triple<DestIterator, DestIterator, DestAccessor> dest)
259 {
260 fourierLogAbsSpectrumInPolarCoordinates(src.first, src.second, src.third, dest.first, dest.second, dest.third);
261 }
262 } //namespace detail
263
264
265
266
267 /********************************************************/
268 /* */
269 /* estimateGlobalRotation */
270 /* */
271 /********************************************************/
272
273 /** \brief Estimate the rotation between two images by means of a normalized cross correlation matching of the FFT spectra.
274
275 This algorithm uses the fast normalized cross correlation to determine a global rotation
276 between two images (from image2 to image1). To derive the rotation, the algorithm performs the following steps:
277 <ol>
278 <li>Transforming both images to the frequency domain using FFTW</li>
279 <li>Create LogAbs PolarCoordinate representations for each spectrum.</li>
280 <li>Determining the final Rotation by performing a fast normalized cross correlation
281 based on the polar representations.</li>
282 </ol>
283 The images are cropped to the corresponding images center-squared before the estimation
284 takes place.
285
286 <b> Declarations:</b>
287
288 <b>\#include</b> \<vigra/affine_registration_fft.hxx\><br>
289 Namespace: vigra
290
291 pass 2D array views:
292 \code
293 namespace vigra {
294 template <class T1, class S1,
295 class T2, class S2>
296 void
297 estimateGlobalRotation(MultiArrayView<2, T1, S1> const & src,
298 MultiArrayView<2, T2, S2> dest,
299 Matrix<double> & affineMatrix,
300 double & correlation_coefficent);
301 }
302 \endcode
303
304 \deprecatedAPI{estimateGlobalRotation}
305 pass \ref ImageIterators and \ref DataAccessors :
306 \code
307 namespace vigra {
308 template <class SrcIterator, class SrcAccessor,
309 class DestIterator, class DestAccessor>
310 void
311 estimateGlobalRotation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
312 DestIterator dul, DestIterator dlr, DestAccessor dest,
313 Matrix<double> & affineMatrix,
314 double & correlation_coefficent)
315 }
316 \endcode
317 use argument objects in conjunction with \ref ArgumentObjectFactories :
318 \code
319 namespace vigra {
320 template <class SrcIterator, class SrcAccessor,
321 class DestIterator, class DestAccessor>
322 void
323 estimateGlobalRotation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
324 triple<DestIterator, DestIterator, DestAccessor> dest,
325 Matrix<double> & affineMatrix,
326 double & correlation_coefficent)
327 }
328 \endcode
329 \deprecatedEnd
330 */
doxygen_overloaded_function(template<...> void estimateGlobalRotation)331 doxygen_overloaded_function(template <...> void estimateGlobalRotation)
332
333 template <class SrcIterator, class SrcAccessor,
334 class DestIterator, class DestAccessor>
335 void
336 estimateGlobalRotation(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
337 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
338 Matrix<double> & affineMatrix,
339 double & correlation_coefficient)
340 {
341 //determine squared centers of both images without consuming any additional memory!
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;
345
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;
349
350 //Determine Shape for united polar coordinate representation
351 Diff2D mean_shape = (s_shape + d_shape)/2;
352
353 int size = min(mean_shape.x, mean_shape.y);
354 if(size %2 == 0)
355 size++;
356
357 const int pc_w = size,
358 pc_h = size*6+1;
359
360
361 DImage s_polCoords(pc_w, pc_h/2),
362 d_polCoords(pc_w, pc_h),
363 ncc(pc_w, pc_h);
364
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));
369
370 detail::fourierLogAbsSpectrumInPolarCoordinates(srcIterRange(d_ul+d2_offset, d_ul+d2_offset+d2_shape, d_acc),
371 destImageRange(d_polCoords));
372
373 //Basic Cross correlation is assumed to be faster here, as there are only pc_h comparisons possible...
374 normalizedCrossCorrelation(srcImageRange(d_polCoords), srcImageRange(s_polCoords), destImage(ncc));
375
376 int max_idx = 0;
377 double max_val = -1;
378
379 const int x=pc_w/2;
380 double val;
381
382 //Only look at a stripe for the maximum angle of rotation
383 //at the image center, at find the best fitting angle...
384 for (int y=0; y<pc_h/2; y++)
385 {
386 val = ncc(x,y+pc_h/4);
387
388 if (val > max_val)
389 {
390 max_idx = y;
391 max_val = val;
392 }
393 }
394
395 double theta = double(max_idx) / pc_h * 360.0;
396
397 affineMatrix = rotationMatrix2DDegrees(theta, vigra::TinyVector<double,2>(s_shape.x/2.0, s_shape.y/2.0));
398 correlation_coefficient = max_val;
399 }
400
401 template <class SrcIterator, class SrcAccessor,
402 class DestIterator, class DestAccessor>
403 inline void
estimateGlobalRotation(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,Matrix<double> & affineMatrix,double & correlation_coefficient)404 estimateGlobalRotation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
405 triple<DestIterator, DestIterator, DestAccessor> dest,
406 Matrix<double> & affineMatrix,
407 double & correlation_coefficient)
408 {
409 estimateGlobalRotation(src.first, src.second, src.third,
410 dest.first, dest.second, dest.third,
411 affineMatrix,
412 correlation_coefficient);
413 }
414
415 template <class T1, class S1,
416 class T2, class S2>
417 inline void
estimateGlobalRotation(MultiArrayView<2,T1,S1> const & src,MultiArrayView<2,T2,S2> dest,Matrix<double> & affineMatrix,double & correlation_coefficent)418 estimateGlobalRotation(MultiArrayView<2, T1, S1> const & src,
419 MultiArrayView<2, T2, S2> dest,
420 Matrix<double> & affineMatrix,
421 double & correlation_coefficent)
422 {
423 estimateGlobalRotation(srcImageRange(src),
424 destImageRange(dest),
425 affineMatrix,
426 correlation_coefficent);
427 }
428
429 /********************************************************/
430 /* */
431 /* estimateGlobalTranslation */
432 /* */
433 /********************************************************/
434
435 /** \brief Estimate the translation between two images by means of a normalized cross correlation matching.
436
437 This algorithm uses the fast normalized cross correlation to determine a global translation
438 between two images (from image2 to image1). To derive the translation, the algorithm consists of differents steps:
439 <ol>
440 <li>Separation of the second image<br/>
441 The second image (the one, for which the translation shall be determined) is cut into five
442 subregions: UpperLeft, UpperRight, Center, LowerLeft and LowerRight, each of 1/4 the size of
443 the original image. Using a border > 0 results in (all) overlapping regions.</li>
444 <li>Cross-Correlation of the subimages to the first image<br/>
445 The subimages are normalized cross-correlated to the (complete) first image.
446 The resulting maximum-likelihood translations and the correlation coefficients are stored for the next step.</li>
447 <li>Determining the final Translation by voting<br/>
448 Each correlation vector gets one vote at the beginning. For each equality of derived motion vectors,
449 the voting to these vectors is incremented. If the maximum number of votes is larger than 1, the vector with the
450 most votes is chosen. If the maximum number of votes is 1, the vector with the maximum likelihood is choosen.</li>
451 </ol>
452 <b> Declarations:</b>
453
454 <b>\#include</b> \<vigra/affine_registration_fft.hxx\><br>
455 Namespace: vigra
456
457 pass 2D array views:
458 \code
459 namespace vigra {
460 template <class T1, class S1,
461 class T2, class S2>
462 void
463 estimateGlobalTranslation(MultiArrayView<2, T1, S1> const & src,
464 MultiArrayView<2, T2, S2> dest,
465 Matrix<double> & affineMatrix,
466 double & correlation_coefficent,
467 Diff2D border = Diff2D(0,0));
468 }
469 \endcode
470
471 \deprecatedAPI{estimateGlobalTranslation}
472 pass \ref ImageIterators and \ref DataAccessors :
473 \code
474 namespace vigra {
475 template <class SrcIterator, class SrcAccessor,
476 class DestIterator, class DestAccessor>
477 void
478 estimateGlobalTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
479 DestIterator dul, DestIterator dlr, DestAccessor dest,
480 Matrix<double> & affineMatrix,
481 double & correlation_coefficent,
482 Diff2D border = Diff2D(0,0))
483 }
484 \endcode
485 use argument objects in conjunction with \ref ArgumentObjectFactories :
486 \code
487 namespace vigra {
488 template <class SrcIterator, class SrcAccessor,
489 class DestIterator, class DestAccessor>
490 void
491 estimateGlobalTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
492 triple<DestIterator, DestIterator, DestAccessor> dest,
493 Matrix<double> & affineMatrix,
494 double & correlation_coefficent,
495 Diff2D border = Diff2D(0,0))
496 }
497 \endcode
498 \deprecatedEnd
499 */
doxygen_overloaded_function(template<...> void estimateGlobalTranslation)500 doxygen_overloaded_function(template <...> void estimateGlobalTranslation)
501
502 template <class SrcIterator, class SrcAccessor,
503 class DestIterator, class DestAccessor>
504 void estimateGlobalTranslation(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
505 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
506 Matrix<double> & affine_matrix,
507 double & correlation_coefficent,
508 Diff2D border = Diff2D(0,0))
509 {
510 ignore_argument(d_lr);
511 typename SrcIterator::difference_type s_shape = s_lr - s_ul;
512
513 //determine matrix by using 5 quater-matches and a maximum likelihood decision:
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--;
517
518 Diff2D q_offsets[5];
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;
524
525 TinyVector<int,2> translation_vectors[5];
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};
528
529 int max_corr_idx=0;
530
531 for (int q=0; q!=5; q++)
532 {
533 Diff2D offset = q_offsets[q];
534
535 //we are searching a transformation from img2 -> transformed image1, thus we switch dest and tmp
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]);
540
541 translation_vectors[q] = translation_vectors[q] - TinyVector<int,2>(offset);
542
543 if(translation_correlations[q] > translation_correlations[max_corr_idx])
544 {
545 max_corr_idx=q;
546 }
547
548 for (int q_old=0; q_old!=q; q_old++)
549 {
550 if (translation_vectors[q] == translation_vectors[q_old])
551 {
552 translation_votes[q_old]++;
553 }
554 }
555 }
556
557 int max_votes_idx=0;
558
559 for (int q=0; q!=5; q++)
560 {
561 if(translation_votes[q] > translation_votes[max_votes_idx])
562 {
563 max_votes_idx=q;
564 }
565 }
566
567 int best_idx=0;
568 if(translation_votes[max_votes_idx] > 1)
569 {
570 best_idx = max_votes_idx;
571 }
572 else
573 {
574 best_idx = max_corr_idx;
575 }
576
577 affine_matrix = translationMatrix2D(translation_vectors[best_idx]);
578 correlation_coefficent = translation_correlations[best_idx];
579 }
580
581 template <class SrcIterator, class SrcAccessor,
582 class DestIterator, class DestAccessor>
583 inline void
estimateGlobalTranslation(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,Matrix<double> & affineMatrix,double & correlation_coefficent,Diff2D border=Diff2D (0,0))584 estimateGlobalTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
585 triple<DestIterator, DestIterator, DestAccessor> dest,
586 Matrix<double> & affineMatrix,
587 double & correlation_coefficent,
588 Diff2D border = Diff2D(0,0))
589 {
590 estimateGlobalTranslation(src.first, src.second, src.third,
591 dest.first, dest.second, dest.third,
592 affineMatrix,
593 correlation_coefficent,
594 border);
595 }
596
597 template <class T1, class S1,
598 class T2, class S2>
599 inline void
estimateGlobalTranslation(MultiArrayView<2,T1,S1> const & src,MultiArrayView<2,T2,S2> dest,Matrix<double> & affineMatrix,double & correlation_coefficent,Diff2D border=Diff2D (0,0))600 estimateGlobalTranslation(MultiArrayView<2, T1, S1> const & src,
601 MultiArrayView<2, T2, S2> dest,
602 Matrix<double> & affineMatrix,
603 double & correlation_coefficent,
604 Diff2D border = Diff2D(0,0))
605 {
606 estimateGlobalTranslation(srcImageRange(src),
607 destImageRange(dest),
608 affineMatrix,
609 correlation_coefficent,
610 border);
611 }
612
613 /********************************************************/
614 /* */
615 /* estimateGlobalRotationTranslation */
616 /* */
617 /********************************************************/
618
619 /** \brief Estimate the (global) rotation and translation between two images by means a normalized cross correlation matching.
620
621 This algorithm use the functions \ref estimateGlobalRotation() and
622 \ref estimateGlobalTranslation() to estimate a matrix which describes
623 the global rotation and translation from the second to the first image.
624
625 <b> Declarations:</b>
626
627 <b>\#include</b> \<vigra/affine_registration_fft.hxx\><br>
628 Namespace: vigra
629
630 pass 2D array views:
631 \code
632 namespace vigra {
633 template <class T1, class S1,
634 class T2, class S2>
635 void
636 estimateGlobalRotationTranslation(MultiArrayView<2, T1, S1> const & src,
637 MultiArrayView<2, T2, S2> dest,
638 Matrix<double> & affineMatrix,
639 double & rotation_correlation,
640 double & translation_correlation,
641 Diff2D border = Diff2D(0,0));
642 }
643 \endcode
644
645 \deprecatedAPI{estimateGlobalRotationTranslation}
646 pass \ref ImageIterators and \ref DataAccessors :
647 \code
648 namespace vigra {
649 template <class SrcIterator, class SrcAccessor,
650 class DestIterator, class DestAccessor>
651 void
652 estimateGlobalRotationTranslation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
653 DestIterator dul, DestIterator dlr, DestAccessor dest,
654 Matrix<double> & affineMatrix,
655 double & rotation_correlation,
656 double & translation_correlation,
657 Diff2D border = Diff2D(0,0))
658 }
659 \endcode
660 use argument objects in conjunction with \ref ArgumentObjectFactories :
661 \code
662 namespace vigra {
663 template <class SrcIterator, class SrcAccessor,
664 class DestIterator, class DestAccessor>
665 void
666 estimateGlobalRotationTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
667 triple<DestIterator, DestIterator, DestAccessor> dest,
668 Matrix<double> & affineMatrix,
669 double & rotation_correlation,
670 double & translation_correlation,
671 Diff2D border = Diff2D(0,0))
672 }
673 \endcode
674 \deprecatedEnd
675 */
doxygen_overloaded_function(template<...> void estimateGlobalRotationTranslation)676 doxygen_overloaded_function(template <...> void estimateGlobalRotationTranslation)
677 template <class SrcIterator, class SrcAccessor,
678 class DestIterator, class DestAccessor>
679 void
680 estimateGlobalRotationTranslation(SrcIterator s_ul, SrcIterator s_lr, SrcAccessor s_acc,
681 DestIterator d_ul, DestIterator d_lr, DestAccessor d_acc,
682 Matrix<double> & affineMatrix,
683 double & rotation_correlation,
684 double & translation_correlation,
685 Diff2D border = Diff2D(0,0))
686 {
687 typename DestIterator::difference_type d_shape = d_lr - d_ul;
688
689 //First step: Estimate rotation from img2 -> img1.
690 Matrix<double> rotation_matrix;
691 estimateGlobalRotation(srcIterRange(s_ul+border, s_lr-border, s_acc),
692 srcIterRange(d_ul+border, d_lr-border, d_acc),
693 rotation_matrix,
694 rotation_correlation);
695
696 //Second step: correct image according to the estimated rotation:
697 FImage tmp(d_shape);
698 SplineImageView<3, double> spl(srcIterRange(s_ul, s_lr, s_acc));
699 affineWarpImage(spl, destImageRange(tmp), rotation_matrix);
700
701 //Third step: find rotation between temp image (of step 2) and dest:
702 Matrix<double> translation_matrix;
703 estimateGlobalTranslation(srcImageRange(tmp),
704 srcIterRange(d_ul, d_lr, d_acc),
705 translation_matrix,
706 translation_correlation,
707 border);
708
709 affineMatrix = rotation_matrix * translation_matrix;
710 }
711
712 template <class SrcIterator, class SrcAccessor,
713 class DestIterator, class DestAccessor>
714 inline void
estimateGlobalRotationTranslation(triple<SrcIterator,SrcIterator,SrcAccessor> src,triple<DestIterator,DestIterator,DestAccessor> dest,Matrix<double> & affineMatrix,double & rotation_correlation,double & translation_correlation,Diff2D border=Diff2D (0,0))715 estimateGlobalRotationTranslation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
716 triple<DestIterator, DestIterator, DestAccessor> dest,
717 Matrix<double> & affineMatrix,
718 double & rotation_correlation,
719 double & translation_correlation,
720 Diff2D border = Diff2D(0,0))
721 {
722 estimateGlobalRotationTranslation(src.first, src.second, src.third,
723 dest.first, dest.second, dest.third,
724 affineMatrix,
725 rotation_correlation,
726 translation_correlation,
727 border);
728 }
729
730 template <class T1, class S1,
731 class T2, class S2>
732 inline void
estimateGlobalRotationTranslation(MultiArrayView<2,T1,S1> const & src,MultiArrayView<2,T2,S2> dest,Matrix<double> & affineMatrix,double & rotation_correlation,double & translation_correlation,Diff2D border=Diff2D (0,0))733 estimateGlobalRotationTranslation(MultiArrayView<2, T1, S1> const & src,
734 MultiArrayView<2, T2, S2> dest,
735 Matrix<double> & affineMatrix,
736 double & rotation_correlation,
737 double & translation_correlation,
738 Diff2D border = Diff2D(0,0))
739 {
740 estimateGlobalRotationTranslation(srcImageRange(src),
741 destImageRange(dest),
742 affineMatrix,
743 rotation_correlation,
744 translation_correlation,
745 border);
746 }
747
748 //@}
749
750 } // namespace vigra
751
752
753 #endif /* VIGRA_AFFINE_REGISTRATION_FFT_HXX */
754