1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2002 by Ullrich Koethe */
4 /* Cognitive Systems Group, University of Hamburg, Germany */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* ( Version 1.5.0, Dec 07 2006 ) */
8 /* The VIGRA Website is */
9 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */
10 /* Please direct questions, bug reports, and contributions to */
11 /* koethe@informatik.uni-hamburg.de or */
12 /* vigra@kogs1.informatik.uni-hamburg.de */
13 /* */
14 /* Permission is hereby granted, free of charge, to any person */
15 /* obtaining a copy of this software and associated documentation */
16 /* files (the "Software"), to deal in the Software without */
17 /* restriction, including without limitation the rights to use, */
18 /* copy, modify, merge, publish, distribute, sublicense, and/or */
19 /* sell copies of the Software, and to permit persons to whom the */
20 /* Software is furnished to do so, subject to the following */
21 /* conditions: */
22 /* */
23 /* The above copyright notice and this permission notice shall be */
24 /* included in all copies or substantial portions of the */
25 /* Software. */
26 /* */
27 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34 /* OTHER DEALINGS IN THE SOFTWARE. */
35 /* */
36 /************************************************************************/
37
38 #ifndef VIGRA_FFTW_HXX
39 #define VIGRA_FFTW_HXX
40
41 #include <cmath>
42 #include <functional>
43 #include "stdimage.hxx"
44 #include "copyimage.hxx"
45 #include "transformimage.hxx"
46 #include "combineimages.hxx"
47 #include "numerictraits.hxx"
48 #include "imagecontainer.hxx"
49 #include <fftw.h>
50
51 namespace vigra {
52
53 /********************************************************/
54 /* */
55 /* FFTWComplex */
56 /* */
57 /********************************************************/
58
59 /* documentation: see fftw3.hxx
60 */
61 class FFTWComplex
62 : public fftw_complex
63 {
64 public:
65 /** The complex' component type, as defined in '<TT>fftw.h</TT>'
66 */
67 typedef fftw_real value_type;
68
69 /** reference type (result of operator[])
70 */
71 typedef fftw_real & reference;
72
73 /** const reference type (result of operator[] const)
74 */
75 typedef fftw_real const & const_reference;
76
77 /** iterator type (result of begin() )
78 */
79 typedef fftw_real * iterator;
80
81 /** const iterator type (result of begin() const)
82 */
83 typedef fftw_real const * const_iterator;
84
85 /** The norm type (result of magnitde())
86 */
87 typedef fftw_real NormType;
88
89 /** The squared norm type (result of squaredMagnitde())
90 */
91 typedef fftw_real SquaredNormType;
92
93 /** Construct from real and imaginary part.
94 Default: 0.
95 */
FFTWComplex(value_type const & ire=0.0,value_type const & iim=0.0)96 FFTWComplex(value_type const & ire = 0.0, value_type const & iim = 0.0)
97 {
98 re = ire;
99 im = iim;
100 }
101
102 /** Copy constructor.
103 */
FFTWComplex(FFTWComplex const & o)104 FFTWComplex(FFTWComplex const & o)
105 : fftw_complex(o)
106 {}
107
108 /** Construct from plain <TT>fftw_complex</TT>.
109 */
FFTWComplex(fftw_complex const & o)110 FFTWComplex(fftw_complex const & o)
111 : fftw_complex(o)
112 {}
113
114 /** Construct from TinyVector.
115 */
116 template <class T>
FFTWComplex(TinyVector<T,2> const & o)117 FFTWComplex(TinyVector<T, 2> const & o)
118 {
119 re = o[0];
120 im = o[1];
121 }
122
123 /** Assignment.
124 */
operator =(FFTWComplex const & o)125 FFTWComplex& operator=(FFTWComplex const & o)
126 {
127 re = o.re;
128 im = o.im;
129 return *this;
130 }
131
132 /** Assignment.
133 */
operator =(fftw_complex const & o)134 FFTWComplex& operator=(fftw_complex const & o)
135 {
136 re = o.re;
137 im = o.im;
138 return *this;
139 }
140
141 /** Assignment.
142 */
operator =(fftw_real const & o)143 FFTWComplex& operator=(fftw_real const & o)
144 {
145 re = o;
146 im = 0.0;
147 return *this;
148 }
149
150 /** Assignment.
151 */
152 template <class T>
operator =(TinyVector<T,2> const & o)153 FFTWComplex& operator=(TinyVector<T, 2> const & o)
154 {
155 re = o[0];
156 im = o[1];
157 return *this;
158 }
159
160 /** Unary negation.
161 */
operator -() const162 FFTWComplex operator-() const
163 { return FFTWComplex(-re, -im); }
164
165 /** Squared magnitude x*conj(x)
166 */
squaredMagnitude() const167 SquaredNormType squaredMagnitude() const
168 { return c_re(*this)*c_re(*this)+c_im(*this)*c_im(*this); }
169
170 /** Magnitude (length of radius vector).
171 */
magnitude() const172 NormType magnitude() const
173 { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
174
175 /** Phase angle.
176 */
phase() const177 value_type phase() const
178 { return VIGRA_CSTD::atan2(c_im(*this),c_re(*this)); }
179
180 /** Access components as if number were a vector.
181 */
operator [](int i)182 reference operator[](int i)
183 { return (&re)[i]; }
184
185 /** Read components as if number were a vector.
186 */
operator [](int i) const187 const_reference operator[](int i) const
188 { return (&re)[i]; }
189
190 /** Length of complex number (always 2).
191 */
size() const192 int size() const
193 { return 2; }
194
begin()195 iterator begin()
196 { return &re; }
197
end()198 iterator end()
199 { return &re + 2; }
200
begin() const201 const_iterator begin() const
202 { return &re; }
203
end() const204 const_iterator end() const
205 { return &re + 2; }
206 };
207
208 /********************************************************/
209 /* */
210 /* FFTWComplex Traits */
211 /* */
212 /********************************************************/
213
214 /* documentation: see fftw3.hxx
215 */
216 template<>
217 struct NumericTraits<fftw_complex>
218 {
219 typedef fftw_complex Type;
220 typedef fftw_complex Promote;
221 typedef fftw_complex RealPromote;
222 typedef fftw_complex ComplexPromote;
223 typedef fftw_real ValueType;
224
225 typedef VigraFalseType isIntegral;
226 typedef VigraFalseType isScalar;
227 typedef NumericTraits<fftw_real>::isSigned isSigned;
228 typedef VigraFalseType isOrdered;
229 typedef VigraTrueType isComplex;
230
zerovigra::NumericTraits231 static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
onevigra::NumericTraits232 static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
nonZerovigra::NumericTraits233 static FFTWComplex nonZero() { return one(); }
234
toPromotevigra::NumericTraits235 static const Promote & toPromote(const Type & v) { return v; }
toRealPromotevigra::NumericTraits236 static const RealPromote & toRealPromote(const Type & v) { return v; }
fromPromotevigra::NumericTraits237 static const Type & fromPromote(const Promote & v) { return v; }
fromRealPromotevigra::NumericTraits238 static const Type & fromRealPromote(const RealPromote & v) { return v; }
239 };
240
241 template<>
242 struct NumericTraits<FFTWComplex>
243 : public NumericTraits<fftw_complex>
244 {
245 typedef FFTWComplex Type;
246 typedef FFTWComplex Promote;
247 typedef FFTWComplex RealPromote;
248 typedef FFTWComplex ComplexPromote;
249 typedef fftw_real ValueType;
250
251 typedef VigraFalseType isIntegral;
252 typedef VigraFalseType isScalar;
253 typedef NumericTraits<fftw_real>::isSigned isSigned;
254 typedef VigraFalseType isOrdered;
255 typedef VigraTrueType isComplex;
256 };
257
258 template<>
259 struct NormTraits<fftw_complex>
260 {
261 typedef fftw_complex Type;
262 typedef fftw_real SquaredNormType;
263 typedef fftw_real NormType;
264 };
265
266 template<>
267 struct NormTraits<FFTWComplex>
268 {
269 typedef FFTWComplex Type;
270 typedef fftw_real SquaredNormType;
271 typedef fftw_real NormType;
272 };
273
274
275 template <>
276 struct PromoteTraits<fftw_complex, fftw_complex>
277 {
278 typedef fftw_complex Promote;
279 };
280
281 template <>
282 struct PromoteTraits<fftw_complex, double>
283 {
284 typedef fftw_complex Promote;
285 };
286
287 template <>
288 struct PromoteTraits<double, fftw_complex>
289 {
290 typedef fftw_complex Promote;
291 };
292
293 template <>
294 struct PromoteTraits<FFTWComplex, FFTWComplex>
295 {
296 typedef FFTWComplex Promote;
297 };
298
299 template <>
300 struct PromoteTraits<FFTWComplex, double>
301 {
302 typedef FFTWComplex Promote;
303 };
304
305 template <>
306 struct PromoteTraits<double, FFTWComplex>
307 {
308 typedef FFTWComplex Promote;
309 };
310
311
312 /********************************************************/
313 /* */
314 /* FFTWComplex Operations */
315 /* */
316 /********************************************************/
317
318 /* documentation: see fftw3.hxx
319 */
operator ==(FFTWComplex const & a,const FFTWComplex & b)320 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) {
321 return c_re(a) == c_re(b) && c_im(a) == c_im(b);
322 }
323
operator !=(FFTWComplex const & a,const FFTWComplex & b)324 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) {
325 return c_re(a) != c_re(b) || c_im(a) != c_im(b);
326 }
327
operator +=(FFTWComplex & a,const FFTWComplex & b)328 inline FFTWComplex &operator +=(FFTWComplex &a, const FFTWComplex &b) {
329 c_re(a) += c_re(b);
330 c_im(a) += c_im(b);
331 return a;
332 }
333
operator -=(FFTWComplex & a,const FFTWComplex & b)334 inline FFTWComplex &operator -=(FFTWComplex &a, const FFTWComplex &b) {
335 c_re(a) -= c_re(b);
336 c_im(a) -= c_im(b);
337 return a;
338 }
339
operator *=(FFTWComplex & a,const FFTWComplex & b)340 inline FFTWComplex &operator *=(FFTWComplex &a, const FFTWComplex &b) {
341 FFTWComplex::value_type t = c_re(a)*c_re(b)-c_im(a)*c_im(b);
342 c_im(a) = c_re(a)*c_im(b)+c_im(a)*c_re(b);
343 c_re(a) = t;
344 return a;
345 }
346
operator /=(FFTWComplex & a,const FFTWComplex & b)347 inline FFTWComplex &operator /=(FFTWComplex &a, const FFTWComplex &b) {
348 FFTWComplex::value_type sm = b.squaredMagnitude();
349 FFTWComplex::value_type t = (c_re(a)*c_re(b)+c_im(a)*c_im(b))/sm;
350 c_im(a) = (c_re(b)*c_im(a)-c_re(a)*c_im(b))/sm;
351 c_re(a) = t;
352 return a;
353 }
354
operator *=(FFTWComplex & a,const double & b)355 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) {
356 c_re(a) *= b;
357 c_im(a) *= b;
358 return a;
359 }
360
operator /=(FFTWComplex & a,const double & b)361 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) {
362 c_re(a) /= b;
363 c_im(a) /= b;
364 return a;
365 }
366
operator +(FFTWComplex a,const FFTWComplex & b)367 inline FFTWComplex operator +(FFTWComplex a, const FFTWComplex &b) {
368 a += b;
369 return a;
370 }
371
operator -(FFTWComplex a,const FFTWComplex & b)372 inline FFTWComplex operator -(FFTWComplex a, const FFTWComplex &b) {
373 a -= b;
374 return a;
375 }
376
operator *(FFTWComplex a,const FFTWComplex & b)377 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) {
378 a *= b;
379 return a;
380 }
381
operator *(FFTWComplex a,const double & b)382 inline FFTWComplex operator *(FFTWComplex a, const double &b) {
383 a *= b;
384 return a;
385 }
386
operator *(const double & a,FFTWComplex b)387 inline FFTWComplex operator *(const double &a, FFTWComplex b) {
388 b *= a;
389 return b;
390 }
391
operator /(FFTWComplex a,const FFTWComplex & b)392 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) {
393 a /= b;
394 return a;
395 }
396
operator /(FFTWComplex a,const double & b)397 inline FFTWComplex operator /(FFTWComplex a, const double &b) {
398 a /= b;
399 return a;
400 }
401
402 using VIGRA_CSTD::abs;
403
abs(const FFTWComplex & a)404 inline FFTWComplex::value_type abs(const FFTWComplex &a)
405 {
406 return a.magnitude();
407 }
408
conj(const FFTWComplex & a)409 inline FFTWComplex conj(const FFTWComplex &a)
410 {
411 return FFTWComplex(a.re, -a.im);
412 }
413
norm(const FFTWComplex & a)414 inline FFTWComplex::NormType norm(const FFTWComplex &a)
415 {
416 return a.magnitude();
417 }
418
squaredNorm(const FFTWComplex & a)419 inline FFTWComplex::SquaredNormType squaredNorm(const FFTWComplex &a)
420 {
421 return a.squaredMagnitude();
422 }
423
424 /********************************************************/
425 /* */
426 /* FFTWRealImage */
427 /* */
428 /********************************************************/
429
430 /* documentation: see fftw3.hxx
431 */
432 typedef BasicImage<fftw_real> FFTWRealImage;
433
434 /********************************************************/
435 /* */
436 /* FFTWComplexImage */
437 /* */
438 /********************************************************/
439
440 template<>
441 struct IteratorTraits<
442 BasicImageIterator<FFTWComplex, FFTWComplex **> >
443 {
444 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> Iterator;
445 typedef Iterator iterator;
446 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> mutable_iterator;
447 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> const_iterator;
448 typedef iterator::iterator_category iterator_category;
449 typedef iterator::value_type value_type;
450 typedef iterator::reference reference;
451 typedef iterator::index_reference index_reference;
452 typedef iterator::pointer pointer;
453 typedef iterator::difference_type difference_type;
454 typedef iterator::row_iterator row_iterator;
455 typedef iterator::column_iterator column_iterator;
456 typedef VectorAccessor<FFTWComplex> default_accessor;
457 typedef VectorAccessor<FFTWComplex> DefaultAccessor;
458 typedef VigraTrueType hasConstantStrides;
459 };
460
461 template<>
462 struct IteratorTraits<
463 ConstBasicImageIterator<FFTWComplex, FFTWComplex **> >
464 {
465 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> Iterator;
466 typedef Iterator iterator;
467 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> mutable_iterator;
468 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> const_iterator;
469 typedef iterator::iterator_category iterator_category;
470 typedef iterator::value_type value_type;
471 typedef iterator::reference reference;
472 typedef iterator::index_reference index_reference;
473 typedef iterator::pointer pointer;
474 typedef iterator::difference_type difference_type;
475 typedef iterator::row_iterator row_iterator;
476 typedef iterator::column_iterator column_iterator;
477 typedef VectorAccessor<FFTWComplex> default_accessor;
478 typedef VectorAccessor<FFTWComplex> DefaultAccessor;
479 typedef VigraTrueType hasConstantStrides;
480 };
481
482 /* documentation: see fftw3.hxx
483 */
484 typedef BasicImage<FFTWComplex> FFTWComplexImage;
485
486 /********************************************************/
487 /* */
488 /* FFTWComplex-Accessors */
489 /* */
490 /********************************************************/
491
492 /* documentation: see fftw3.hxx
493 */
494 class FFTWRealAccessor
495 {
496 public:
497
498 /// The accessor's value type.
499 typedef fftw_real value_type;
500
501 /// Read real part at iterator position.
502 template <class ITERATOR>
operator ()(ITERATOR const & i) const503 value_type operator()(ITERATOR const & i) const {
504 return c_re(*i);
505 }
506
507 /// Read real part at offset from iterator position.
508 template <class ITERATOR, class DIFFERENCE>
operator ()(ITERATOR const & i,DIFFERENCE d) const509 value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
510 return c_re(i[d]);
511 }
512
513 /// Write real part at iterator position.
514 template <class ITERATOR>
set(value_type const & v,ITERATOR const & i) const515 void set(value_type const & v, ITERATOR const & i) const {
516 c_re(*i)= v;
517 }
518
519 /// Write real part at offset from iterator position.
520 template <class ITERATOR, class DIFFERENCE>
set(value_type const & v,ITERATOR const & i,DIFFERENCE d) const521 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
522 c_re(i[d])= v;
523 }
524 };
525
526 /* documentation: see fftw3.hxx
527 */
528 class FFTWImaginaryAccessor
529 {
530 public:
531 /// The accessor's value type.
532 typedef fftw_real value_type;
533
534 /// Read imaginary part at iterator position.
535 template <class ITERATOR>
operator ()(ITERATOR const & i) const536 value_type operator()(ITERATOR const & i) const {
537 return c_im(*i);
538 }
539
540 /// Read imaginary part at offset from iterator position.
541 template <class ITERATOR, class DIFFERENCE>
operator ()(ITERATOR const & i,DIFFERENCE d) const542 value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
543 return c_im(i[d]);
544 }
545
546 /// Write imaginary part at iterator position.
547 template <class ITERATOR>
set(value_type const & v,ITERATOR const & i) const548 void set(value_type const & v, ITERATOR const & i) const {
549 c_im(*i)= v;
550 }
551
552 /// Write imaginary part at offset from iterator position.
553 template <class ITERATOR, class DIFFERENCE>
set(value_type const & v,ITERATOR const & i,DIFFERENCE d) const554 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
555 c_im(i[d])= v;
556 }
557 };
558
559 /* documentation: see fftw3.hxx
560 */
561 class FFTWWriteRealAccessor: public FFTWRealAccessor
562 {
563 public:
564 /// The accessor's value type.
565 typedef fftw_real value_type;
566
567 /** Write real number at iterator position. Set imaginary part
568 to 0.
569 */
570 template <class ITERATOR>
set(value_type const & v,ITERATOR const & i) const571 void set(value_type const & v, ITERATOR const & i) const {
572 c_re(*i)= v;
573 c_im(*i)= 0;
574 }
575
576 /** Write real number at offset from iterator position. Set imaginary part
577 to 0.
578 */
579 template <class ITERATOR, class DIFFERENCE>
set(value_type const & v,ITERATOR const & i,DIFFERENCE d) const580 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
581 c_re(i[d])= v;
582 c_im(i[d])= 0;
583 }
584 };
585
586 /* documentation: see fftw3.hxx
587 */
588 class FFTWMagnitudeAccessor
589 {
590 public:
591 /// The accessor's value type.
592 typedef fftw_real value_type;
593
594 /// Read magnitude at iterator position.
595 template <class ITERATOR>
operator ()(ITERATOR const & i) const596 value_type operator()(ITERATOR const & i) const {
597 return (*i).magnitude();
598 }
599
600 /// Read magnitude at offset from iterator position.
601 template <class ITERATOR, class DIFFERENCE>
operator ()(ITERATOR const & i,DIFFERENCE d) const602 value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
603 return (i[d]).magnitude();
604 }
605 };
606
607 /* documentation: see fftw3.hxx
608 */
609 class FFTWPhaseAccessor
610 {
611 public:
612 /// The accessor's value type.
613 typedef fftw_real value_type;
614
615 /// Read phase at iterator position.
616 template <class ITERATOR>
operator ()(ITERATOR const & i) const617 value_type operator()(ITERATOR const & i) const {
618 return (*i).phase();
619 }
620
621 /// Read phase at offset from iterator position.
622 template <class ITERATOR, class DIFFERENCE>
operator ()(ITERATOR const & i,DIFFERENCE d) const623 value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
624 return (i[d]).phase();
625 }
626 };
627
628 /********************************************************/
629 /* */
630 /* Fourier Transform */
631 /* */
632 /********************************************************/
633
634 /** \page FourierTransformFFTW2 Fast Fourier Transform
635
636 This documentation describes the deprecated VIGRA interface to
637 FFTW 2. Use the \link FourierTransform interface to the newer
638 version FFTW 3\endlink instead.
639
640 VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
641 Transform</a> package to perform Fourier transformations. VIGRA
642 provides a wrapper for FFTW's complex number type (FFTWComplex),
643 but FFTW's functions are used verbatim. If the image is stored as
644 a FFTWComplexImage, a FFT is performed like this:
645
646 \code
647 vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
648 ... // fill image with data
649
650 // create a plan for optimal performance
651 fftwnd_plan forwardPlan=
652 fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE );
653
654 // calculate FFT
655 fftwnd_one(forwardPlan, spatial.begin(), fourier.begin());
656 \endcode
657
658 Note that in the creation of a plan, the height must be given
659 first. Note also that <TT>spatial.begin()</TT> may only be passed
660 to <TT>fftwnd_one</TT> if the transform shall be applied to the
661 entire image. When you want to retrict operation to an ROI, you
662 create a copy of the ROI in an image of appropriate size.
663
664 More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
665
666 FFTW produces fourier images that have the DC component (the
667 origin of the Fourier space) in the upper left corner. Often, one
668 wants the origin in the center of the image, so that frequencies
669 always increase towards the border of the image. This can be
670 achieved by calling \ref moveDCToCenter(). The inverse
671 transformation is done by \ref moveDCToUpperLeft().
672
673 <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
674 Namespace: vigra
675 */
676
677 /********************************************************/
678 /* */
679 /* moveDCToCenter */
680 /* */
681 /********************************************************/
682
683 /* documentation: see fftw3.hxx
684 */
685 template <class SrcImageIterator, class SrcAccessor,
686 class DestImageIterator, class DestAccessor>
moveDCToCenter(SrcImageIterator src_upperleft,SrcImageIterator src_lowerright,SrcAccessor sa,DestImageIterator dest_upperleft,DestAccessor da)687 void moveDCToCenter(SrcImageIterator src_upperleft,
688 SrcImageIterator src_lowerright, SrcAccessor sa,
689 DestImageIterator dest_upperleft, DestAccessor da)
690 {
691 int w= src_lowerright.x - src_upperleft.x;
692 int h= src_lowerright.y - src_upperleft.y;
693 int w1 = w/2;
694 int h1 = h/2;
695 int w2 = (w+1)/2;
696 int h2 = (h+1)/2;
697
698 // 2. Quadrant zum 4.
699 copyImage(srcIterRange(src_upperleft,
700 src_upperleft + Diff2D(w2, h2), sa),
701 destIter (dest_upperleft + Diff2D(w1, h1), da));
702
703 // 4. Quadrant zum 2.
704 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
705 src_lowerright, sa),
706 destIter (dest_upperleft, da));
707
708 // 1. Quadrant zum 3.
709 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
710 src_upperleft + Diff2D(w, h2), sa),
711 destIter (dest_upperleft + Diff2D(0, h1), da));
712
713 // 3. Quadrant zum 1.
714 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
715 src_upperleft + Diff2D(w2, h), sa),
716 destIter (dest_upperleft + Diff2D(w1, 0), da));
717 }
718
719 template <class SrcImageIterator, class SrcAccessor,
720 class DestImageIterator, class DestAccessor>
moveDCToCenter(triple<SrcImageIterator,SrcImageIterator,SrcAccessor> src,pair<DestImageIterator,DestAccessor> dest)721 inline void moveDCToCenter(
722 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
723 pair<DestImageIterator, DestAccessor> dest)
724 {
725 moveDCToCenter(src.first, src.second, src.third,
726 dest.first, dest.second);
727 }
728
729 /********************************************************/
730 /* */
731 /* moveDCToUpperLeft */
732 /* */
733 /********************************************************/
734
735 /* documentation: see fftw3.hxx
736 */
737 template <class SrcImageIterator, class SrcAccessor,
738 class DestImageIterator, class DestAccessor>
moveDCToUpperLeft(SrcImageIterator src_upperleft,SrcImageIterator src_lowerright,SrcAccessor sa,DestImageIterator dest_upperleft,DestAccessor da)739 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
740 SrcImageIterator src_lowerright, SrcAccessor sa,
741 DestImageIterator dest_upperleft, DestAccessor da)
742 {
743 int w= src_lowerright.x - src_upperleft.x;
744 int h= src_lowerright.y - src_upperleft.y;
745 int w2 = w/2;
746 int h2 = h/2;
747 int w1 = (w+1)/2;
748 int h1 = (h+1)/2;
749
750 // 2. Quadrant zum 4.
751 copyImage(srcIterRange(src_upperleft,
752 src_upperleft + Diff2D(w2, h2), sa),
753 destIter (dest_upperleft + Diff2D(w1, h1), da));
754
755 // 4. Quadrant zum 2.
756 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
757 src_lowerright, sa),
758 destIter (dest_upperleft, da));
759
760 // 1. Quadrant zum 3.
761 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
762 src_upperleft + Diff2D(w, h2), sa),
763 destIter (dest_upperleft + Diff2D(0, h1), da));
764
765 // 3. Quadrant zum 1.
766 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
767 src_upperleft + Diff2D(w2, h), sa),
768 destIter (dest_upperleft + Diff2D(w1, 0), da));
769 }
770
771 template <class SrcImageIterator, class SrcAccessor,
772 class DestImageIterator, class DestAccessor>
moveDCToUpperLeft(triple<SrcImageIterator,SrcImageIterator,SrcAccessor> src,pair<DestImageIterator,DestAccessor> dest)773 inline void moveDCToUpperLeft(
774 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
775 pair<DestImageIterator, DestAccessor> dest)
776 {
777 moveDCToUpperLeft(src.first, src.second, src.third,
778 dest.first, dest.second);
779 }
780
781 /********************************************************/
782 /* */
783 /* applyFourierFilter */
784 /* */
785 /********************************************************/
786
787 /* documentation: see fftw3.hxx
788 */
789
790 // applyFourierFilter versions without fftwnd_plans:
791 template <class SrcImageIterator, class SrcAccessor,
792 class FilterImageIterator, class FilterAccessor,
793 class DestImageIterator, class DestAccessor>
794 inline
applyFourierFilter(triple<SrcImageIterator,SrcImageIterator,SrcAccessor> src,pair<FilterImageIterator,FilterAccessor> filter,pair<DestImageIterator,DestAccessor> dest)795 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
796 pair<FilterImageIterator, FilterAccessor> filter,
797 pair<DestImageIterator, DestAccessor> dest)
798 {
799 applyFourierFilter(src.first, src.second, src.third,
800 filter.first, filter.second,
801 dest.first, dest.second);
802 }
803
804 template <class SrcImageIterator, class SrcAccessor,
805 class FilterImageIterator, class FilterAccessor,
806 class DestImageIterator, class DestAccessor>
applyFourierFilter(SrcImageIterator srcUpperLeft,SrcImageIterator srcLowerRight,SrcAccessor sa,FilterImageIterator filterUpperLeft,FilterAccessor fa,DestImageIterator destUpperLeft,DestAccessor da)807 void applyFourierFilter(SrcImageIterator srcUpperLeft,
808 SrcImageIterator srcLowerRight, SrcAccessor sa,
809 FilterImageIterator filterUpperLeft, FilterAccessor fa,
810 DestImageIterator destUpperLeft, DestAccessor da)
811 {
812 // copy real input images into a complex one...
813 int w= srcLowerRight.x - srcUpperLeft.x;
814 int h= srcLowerRight.y - srcUpperLeft.y;
815
816 FFTWComplexImage workImage(w, h);
817 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
818 destImage(workImage, FFTWWriteRealAccessor()));
819
820 // ...and call the impl
821 FFTWComplexImage const & cworkImage = workImage;
822 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
823 filterUpperLeft, fa,
824 destUpperLeft, da);
825 }
826
827 typedef FFTWComplexImage::const_traverser FFTWConstTraverser;
828
829 template <class FilterImageIterator, class FilterAccessor,
830 class DestImageIterator, class DestAccessor>
831 inline
applyFourierFilter(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,FilterImageIterator filterUpperLeft,FilterAccessor fa,DestImageIterator destUpperLeft,DestAccessor da)832 void applyFourierFilter(
833 FFTWComplexImage::const_traverser srcUpperLeft,
834 FFTWComplexImage::const_traverser srcLowerRight,
835 FFTWComplexImage::ConstAccessor sa,
836 FilterImageIterator filterUpperLeft, FilterAccessor fa,
837 DestImageIterator destUpperLeft, DestAccessor da)
838 {
839 int w= srcLowerRight.x - srcUpperLeft.x;
840 int h= srcLowerRight.y - srcUpperLeft.y;
841
842 // test for right memory layout (fftw expects a 2*width*height floats array)
843 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
844 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
845 filterUpperLeft, fa,
846 destUpperLeft, da);
847 else
848 {
849 FFTWComplexImage workImage(w, h);
850 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
851 destImage(workImage));
852
853 FFTWComplexImage const & cworkImage = workImage;
854 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
855 filterUpperLeft, fa,
856 destUpperLeft, da);
857 }
858 }
859
860 template <class FilterImageIterator, class FilterAccessor,
861 class DestImageIterator, class DestAccessor>
applyFourierFilterImpl(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,FilterImageIterator filterUpperLeft,FilterAccessor fa,DestImageIterator destUpperLeft,DestAccessor da)862 void applyFourierFilterImpl(
863 FFTWComplexImage::const_traverser srcUpperLeft,
864 FFTWComplexImage::const_traverser srcLowerRight,
865 FFTWComplexImage::ConstAccessor sa,
866 FilterImageIterator filterUpperLeft, FilterAccessor fa,
867 DestImageIterator destUpperLeft, DestAccessor da)
868 {
869 // create plans and call variant with plan parameters
870 int w= srcLowerRight.x - srcUpperLeft.x;
871 int h= srcLowerRight.y - srcUpperLeft.y;
872
873 fftwnd_plan forwardPlan=
874 fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE );
875 fftwnd_plan backwardPlan=
876 fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
877
878 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
879 filterUpperLeft, fa,
880 destUpperLeft, da,
881 forwardPlan, backwardPlan);
882
883 fftwnd_destroy_plan(forwardPlan);
884 fftwnd_destroy_plan(backwardPlan);
885 }
886
887 // applyFourierFilter versions with fftwnd_plans:
888 template <class SrcImageIterator, class SrcAccessor,
889 class FilterImageIterator, class FilterAccessor,
890 class DestImageIterator, class DestAccessor>
891 inline
applyFourierFilter(triple<SrcImageIterator,SrcImageIterator,SrcAccessor> src,pair<FilterImageIterator,FilterAccessor> filter,pair<DestImageIterator,DestAccessor> dest,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)892 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
893 pair<FilterImageIterator, FilterAccessor> filter,
894 pair<DestImageIterator, DestAccessor> dest,
895 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
896 {
897 applyFourierFilter(src.first, src.second, src.third,
898 filter.first, filter.second,
899 dest.first, dest.second,
900 forwardPlan, backwardPlan);
901 }
902
903 template <class SrcImageIterator, class SrcAccessor,
904 class FilterImageIterator, class FilterAccessor,
905 class DestImageIterator, class DestAccessor>
applyFourierFilter(SrcImageIterator srcUpperLeft,SrcImageIterator srcLowerRight,SrcAccessor sa,FilterImageIterator filterUpperLeft,FilterAccessor fa,DestImageIterator destUpperLeft,DestAccessor da,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)906 void applyFourierFilter(SrcImageIterator srcUpperLeft,
907 SrcImageIterator srcLowerRight, SrcAccessor sa,
908 FilterImageIterator filterUpperLeft, FilterAccessor fa,
909 DestImageIterator destUpperLeft, DestAccessor da,
910 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
911 {
912 int w= srcLowerRight.x - srcUpperLeft.x;
913 int h= srcLowerRight.y - srcUpperLeft.y;
914
915 FFTWComplexImage workImage(w, h);
916 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
917 destImage(workImage, FFTWWriteRealAccessor()));
918
919 FFTWComplexImage const & cworkImage = workImage;
920 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
921 filterUpperLeft, fa,
922 destUpperLeft, da,
923 forwardPlan, backwardPlan);
924 }
925
926 template <class FilterImageIterator, class FilterAccessor,
927 class DestImageIterator, class DestAccessor>
928 inline
applyFourierFilter(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,FilterImageIterator filterUpperLeft,FilterAccessor fa,DestImageIterator destUpperLeft,DestAccessor da,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)929 void applyFourierFilter(
930 FFTWComplexImage::const_traverser srcUpperLeft,
931 FFTWComplexImage::const_traverser srcLowerRight,
932 FFTWComplexImage::ConstAccessor sa,
933 FilterImageIterator filterUpperLeft, FilterAccessor fa,
934 DestImageIterator destUpperLeft, DestAccessor da,
935 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
936 {
937 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
938 filterUpperLeft, fa,
939 destUpperLeft, da,
940 forwardPlan, backwardPlan);
941 }
942
943 template <class FilterImageIterator, class FilterAccessor,
944 class DestImageIterator, class DestAccessor>
applyFourierFilterImpl(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,FilterImageIterator filterUpperLeft,FilterAccessor fa,DestImageIterator destUpperLeft,DestAccessor da,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)945 void applyFourierFilterImpl(
946 FFTWComplexImage::const_traverser srcUpperLeft,
947 FFTWComplexImage::const_traverser srcLowerRight,
948 FFTWComplexImage::ConstAccessor sa,
949 FilterImageIterator filterUpperLeft, FilterAccessor fa,
950 DestImageIterator destUpperLeft, DestAccessor da,
951 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
952 {
953 FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
954
955 // FFT from srcImage to complexResultImg
956 fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)),
957 complexResultImg.begin());
958
959 // convolve in freq. domain (in complexResultImg)
960 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
961 destImage(complexResultImg), std::multiplies<FFTWComplex>());
962
963 // FFT back into spatial domain (inplace in complexResultImg)
964 fftwnd_one(backwardPlan, complexResultImg.begin(), 0);
965
966 typedef typename
967 NumericTraits<typename DestAccessor::value_type>::isScalar
968 isScalarResult;
969
970 // normalization (after FFTs), maybe stripping imaginary part
971 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
972 isScalarResult());
973 }
974
975 template <class DestImageIterator, class DestAccessor>
applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,DestImageIterator destUpperLeft,DestAccessor da,VigraFalseType)976 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
977 DestImageIterator destUpperLeft,
978 DestAccessor da,
979 VigraFalseType)
980 {
981 double normFactor= 1.0/(srcImage.width() * srcImage.height());
982
983 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
984 {
985 DestImageIterator dIt= destUpperLeft;
986 for(int x= 0; x< srcImage.width(); x++, dIt.x++)
987 {
988 da.setComponent(srcImage(x, y).re*normFactor, dIt, 0);
989 da.setComponent(srcImage(x, y).im*normFactor, dIt, 1);
990 }
991 }
992 }
993
994 inline
applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,FFTWComplexImage::traverser destUpperLeft,FFTWComplexImage::Accessor da,VigraFalseType)995 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
996 FFTWComplexImage::traverser destUpperLeft,
997 FFTWComplexImage::Accessor da,
998 VigraFalseType)
999 {
1000 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
1001 linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height())));
1002 }
1003
1004 template <class DestImageIterator, class DestAccessor>
applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,DestImageIterator destUpperLeft,DestAccessor da,VigraTrueType)1005 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
1006 DestImageIterator destUpperLeft,
1007 DestAccessor da,
1008 VigraTrueType)
1009 {
1010 double normFactor= 1.0/(srcImage.width() * srcImage.height());
1011
1012 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
1013 {
1014 DestImageIterator dIt= destUpperLeft;
1015 for(int x= 0; x< srcImage.width(); x++, dIt.x++)
1016 da.set(srcImage(x, y).re*normFactor, dIt);
1017 }
1018 }
1019
1020 /**********************************************************/
1021 /* */
1022 /* applyFourierFilterFamily */
1023 /* */
1024 /**********************************************************/
1025
1026 /* documentation: see fftw3.hxx
1027 */
1028
1029 // applyFourierFilterFamily versions without fftwnd_plans:
1030 template <class SrcImageIterator, class SrcAccessor,
1031 class FilterType, class DestImage>
1032 inline
applyFourierFilterFamily(triple<SrcImageIterator,SrcImageIterator,SrcAccessor> src,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results)1033 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1034 const ImageArray<FilterType> &filters,
1035 ImageArray<DestImage> &results)
1036 {
1037 applyFourierFilterFamily(src.first, src.second, src.third,
1038 filters, results);
1039 }
1040
1041 template <class SrcImageIterator, class SrcAccessor,
1042 class FilterType, class DestImage>
applyFourierFilterFamily(SrcImageIterator srcUpperLeft,SrcImageIterator srcLowerRight,SrcAccessor sa,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results)1043 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
1044 SrcImageIterator srcLowerRight, SrcAccessor sa,
1045 const ImageArray<FilterType> &filters,
1046 ImageArray<DestImage> &results)
1047 {
1048 int w= srcLowerRight.x - srcUpperLeft.x;
1049 int h= srcLowerRight.y - srcUpperLeft.y;
1050
1051 FFTWComplexImage workImage(w, h);
1052 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1053 destImage(workImage, FFTWWriteRealAccessor()));
1054
1055 FFTWComplexImage const & cworkImage = workImage;
1056 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1057 filters, results);
1058 }
1059
1060 template <class FilterType, class DestImage>
1061 inline
applyFourierFilterFamily(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results)1062 void applyFourierFilterFamily(
1063 FFTWComplexImage::const_traverser srcUpperLeft,
1064 FFTWComplexImage::const_traverser srcLowerRight,
1065 FFTWComplexImage::ConstAccessor sa,
1066 const ImageArray<FilterType> &filters,
1067 ImageArray<DestImage> &results)
1068 {
1069 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
1070 filters, results);
1071 }
1072
1073 template <class FilterType, class DestImage>
applyFourierFilterFamilyImpl(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results)1074 void applyFourierFilterFamilyImpl(
1075 FFTWComplexImage::const_traverser srcUpperLeft,
1076 FFTWComplexImage::const_traverser srcLowerRight,
1077 FFTWComplexImage::ConstAccessor sa,
1078 const ImageArray<FilterType> &filters,
1079 ImageArray<DestImage> &results)
1080 {
1081 int w= srcLowerRight.x - srcUpperLeft.x;
1082 int h= srcLowerRight.y - srcUpperLeft.y;
1083
1084 fftwnd_plan forwardPlan=
1085 fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE );
1086 fftwnd_plan backwardPlan=
1087 fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
1088
1089 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
1090 filters, results,
1091 forwardPlan, backwardPlan);
1092
1093 fftwnd_destroy_plan(forwardPlan);
1094 fftwnd_destroy_plan(backwardPlan);
1095 }
1096
1097 // applyFourierFilterFamily versions with fftwnd_plans:
1098 template <class SrcImageIterator, class SrcAccessor,
1099 class FilterType, class DestImage>
1100 inline
applyFourierFilterFamily(triple<SrcImageIterator,SrcImageIterator,SrcAccessor> src,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)1101 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1102 const ImageArray<FilterType> &filters,
1103 ImageArray<DestImage> &results,
1104 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
1105 {
1106 applyFourierFilterFamily(src.first, src.second, src.third,
1107 filters, results,
1108 forwardPlan, backwardPlan);
1109 }
1110
1111 template <class SrcImageIterator, class SrcAccessor,
1112 class FilterType, class DestImage>
applyFourierFilterFamily(SrcImageIterator srcUpperLeft,SrcImageIterator srcLowerRight,SrcAccessor sa,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)1113 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
1114 SrcImageIterator srcLowerRight, SrcAccessor sa,
1115 const ImageArray<FilterType> &filters,
1116 ImageArray<DestImage> &results,
1117 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
1118 {
1119 int w= srcLowerRight.x - srcUpperLeft.x;
1120 int h= srcLowerRight.y - srcUpperLeft.y;
1121
1122 FFTWComplexImage workImage(w, h);
1123 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1124 destImage(workImage, FFTWWriteRealAccessor()));
1125
1126 FFTWComplexImage const & cworkImage = workImage;
1127 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1128 filters, results,
1129 forwardPlan, backwardPlan);
1130 }
1131
1132 template <class FilterType, class DestImage>
1133 inline
applyFourierFilterFamily(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)1134 void applyFourierFilterFamily(
1135 FFTWComplexImage::const_traverser srcUpperLeft,
1136 FFTWComplexImage::const_traverser srcLowerRight,
1137 FFTWComplexImage::ConstAccessor sa,
1138 const ImageArray<FilterType> &filters,
1139 ImageArray<DestImage> &results,
1140 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
1141 {
1142 int w= srcLowerRight.x - srcUpperLeft.x;
1143 int h= srcLowerRight.y - srcUpperLeft.y;
1144
1145 // test for right memory layout (fftw expects a 2*width*height floats array)
1146 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
1147 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
1148 filters, results,
1149 forwardPlan, backwardPlan);
1150 else
1151 {
1152 FFTWComplexImage workImage(w, h);
1153 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1154 destImage(workImage));
1155
1156 FFTWComplexImage const & cworkImage = workImage;
1157 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1158 filters, results,
1159 forwardPlan, backwardPlan);
1160 }
1161 }
1162
1163 template <class FilterType, class DestImage>
applyFourierFilterFamilyImpl(FFTWComplexImage::const_traverser srcUpperLeft,FFTWComplexImage::const_traverser srcLowerRight,FFTWComplexImage::ConstAccessor sa,const ImageArray<FilterType> & filters,ImageArray<DestImage> & results,const fftwnd_plan & forwardPlan,const fftwnd_plan & backwardPlan)1164 void applyFourierFilterFamilyImpl(
1165 FFTWComplexImage::const_traverser srcUpperLeft,
1166 FFTWComplexImage::const_traverser srcLowerRight,
1167 FFTWComplexImage::ConstAccessor sa,
1168 const ImageArray<FilterType> &filters,
1169 ImageArray<DestImage> &results,
1170 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
1171 {
1172 // make sure the filter images have the right dimensions
1173 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
1174 "applyFourierFilterFamily called with src image size != filters.imageSize()!");
1175
1176 // make sure the result image array has the right dimensions
1177 results.resize(filters.size());
1178 results.resizeImages(filters.imageSize());
1179
1180 // FFT from srcImage to freqImage
1181 int w= srcLowerRight.x - srcUpperLeft.x;
1182 int h= srcLowerRight.y - srcUpperLeft.y;
1183
1184 FFTWComplexImage freqImage(w, h);
1185 FFTWComplexImage result(w, h);
1186
1187 fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)), freqImage.begin());
1188
1189 typedef typename
1190 NumericTraits<typename DestImage::Accessor::value_type>::isScalar
1191 isScalarResult;
1192
1193 // convolve with filters in freq. domain
1194 for (unsigned int i= 0; i < filters.size(); i++)
1195 {
1196 combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
1197 destImage(result), std::multiplies<FFTWComplex>());
1198
1199 // FFT back into spatial domain (inplace in destImage)
1200 fftwnd_one(backwardPlan, result.begin(), 0);
1201
1202 // normalization (after FFTs), maybe stripping imaginary part
1203 applyFourierFilterImplNormalization(result,
1204 results[i].upperLeft(), results[i].accessor(),
1205 isScalarResult());
1206 }
1207 }
1208
1209 } // namespace vigra
1210
1211 #endif // VIGRA_FFTW_HXX
1212