[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

fftw3.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
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_FFTW3_HXX
37 #define VIGRA_FFTW3_HXX
38 
39 #include <cmath>
40 #include <functional>
41 #include "stdimage.hxx"
42 #include "copyimage.hxx"
43 #include "transformimage.hxx"
44 #include "combineimages.hxx"
45 #include "numerictraits.hxx"
46 #include "imagecontainer.hxx"
47 #include <fftw3.h>
48 
49 namespace vigra {
50 
51 typedef double fftw_real;
52 
53 /********************************************************/
54 /* */
55 /* FFTWComplex */
56 /* */
57 /********************************************************/
58 
59 /** \brief Wrapper class for the FFTW type '<TT>fftw_complex</TT>'.
60 
61  This class provides constructors and other member functions
62  for the C struct '<TT>fftw_complex</TT>'. This struct is the basic
63  pixel type of the <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a>
64  library. It inherits the data members '<TT>re</TT>' and '<TT>im</TT>'
65  that denote the real and imaginary part of the number. In addition it
66  defines transformations to polar coordinates,
67  as well as \ref FFTWComplexOperators "arithmetic operators" and
68  \ref FFTWComplexAccessors "accessors".
69 
70  FFTWComplex implements the concepts \ref AlgebraicField and
71  \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt>
72  and <tt>FFTWComplexImage</tt> are defined.
73 
74  <b>See also:</b>
75  <ul>
76  <li> \ref FFTWComplexTraits
77  <li> \ref FFTWComplexOperators
78  <li> \ref FFTWComplexAccessors
79  </ul>
80 
81  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
82  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
83  Namespace: vigra
84 */
86 {
87  fftw_complex data_;
88 
89  public:
90  /** The complex' component type, as defined in '<TT>fftw3.h</TT>'
91  */
92  typedef fftw_real value_type;
93 
94  /** reference type (result of operator[])
95  */
96  typedef fftw_real & reference;
97 
98  /** const reference type (result of operator[] const)
99  */
100  typedef fftw_real const & const_reference;
101 
102  /** iterator type (result of begin() )
103  */
104  typedef fftw_real * iterator;
105 
106  /** const iterator type (result of begin() const)
107  */
108  typedef fftw_real const * const_iterator;
109 
110  /** The norm type (result of magnitde())
111  */
112  typedef fftw_real NormType;
113 
114  /** The squared norm type (result of squaredMagnitde())
115  */
116  typedef fftw_real SquaredNormType;
117 
118  /** Construct from real and imaginary part.
119  Default: 0.
120  */
121  FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0)
122  {
123  data_[0] = re;
124  data_[1] = im;
125  }
126 
127  /** Copy constructor.
128  */
130  {
131  data_[0] = o.data_[0];
132  data_[1] = o.data_[1];
133  }
134 
135  /** Construct from plain <TT>fftw_complex</TT>.
136  */
137  FFTWComplex(fftw_complex const & o)
138  {
139  data_[0] = o[0];
140  data_[1] = o[1];
141  }
142 
143  /** Construct from TinyVector.
144  */
145  template <class T>
147  {
148  data_[0] = o[0];
149  data_[1] = o[1];
150  }
151 
152  /** Assignment.
153  */
155  {
156  data_[0] = o.data_[0];
157  data_[1] = o.data_[1];
158  return *this;
159  }
160 
161  /** Assignment.
162  */
163  FFTWComplex& operator=(fftw_complex const & o)
164  {
165  data_[0] = o[0];
166  data_[1] = o[1];
167  return *this;
168  }
169 
170  /** Assignment.
171  */
172  FFTWComplex& operator=(fftw_real const & o)
173  {
174  data_[0] = o;
175  data_[1] = 0.0;
176  return *this;
177  }
178 
179  /** Assignment.
180  */
181  template <class T>
183  {
184  data_[0] = o[0];
185  data_[1] = o[1];
186  return *this;
187  }
188 
189  reference re()
190  { return data_[0]; }
191 
192  const_reference re() const
193  { return data_[0]; }
194 
195  reference im()
196  { return data_[1]; }
197 
198  const_reference im() const
199  { return data_[1]; }
200 
201  /** Unary negation.
202  */
204  { return FFTWComplex(-data_[0], -data_[1]); }
205 
206  /** Squared magnitude x*conj(x)
207  */
209  { return data_[0]*data_[0]+data_[1]*data_[1]; }
210 
211  /** Magnitude (length of radius vector).
212  */
214  { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
215 
216  /** Phase angle.
217  */
219  { return VIGRA_CSTD::atan2(data_[1], data_[0]); }
220 
221  /** Access components as if number were a vector.
222  */
224  { return data_[i]; }
225 
226  /** Read components as if number were a vector.
227  */
229  { return data_[i]; }
230 
231  /** Length of complex number (always 2).
232  */
233  int size() const
234  { return 2; }
235 
236  iterator begin()
237  { return data_; }
238 
239  iterator end()
240  { return data_ + 2; }
241 
242  const_iterator begin() const
243  { return data_; }
244 
245  const_iterator end() const
246  { return data_ + 2; }
247 };
248 
249 /********************************************************/
250 /* */
251 /* FFTWComplexTraits */
252 /* */
253 /********************************************************/
254 
255 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex
256 
257  The numeric and promote traits for fftw_complex and FFTWComplex follow
258  the general specifications for \ref NumericPromotionTraits and
259  \ref AlgebraicField. They are explicitly specialized for the types
260  involved:
261 
262  \code
263 
264  template<>
265  struct NumericTraits<fftw_complex>
266  {
267  typedef fftw_complex Promote;
268  typedef fftw_complex RealPromote;
269  typedef fftw_complex ComplexPromote;
270  typedef fftw_real ValueType;
271 
272  typedef VigraFalseType isIntegral;
273  typedef VigraFalseType isScalar;
274  typedef VigraFalseType isOrdered;
275  typedef VigraTrueType isComplex;
276 
277  // etc.
278  };
279 
280  template<>
281  struct NumericTraits<FFTWComplex>
282  {
283  typedef FFTWComplex Promote;
284  typedef FFTWComplex RealPromote;
285  typedef FFTWComplex ComplexPromote;
286  typedef fftw_real ValueType;
287 
288  typedef VigraFalseType isIntegral;
289  typedef VigraFalseType isScalar;
290  typedef VigraFalseType isOrdered;
291  typedef VigraTrueType isComplex;
292 
293  // etc.
294  };
295 
296  template<>
297  struct NormTraits<fftw_complex>
298  {
299  typedef fftw_complex Type;
300  typedef fftw_real SquaredNormType;
301  typedef fftw_real NormType;
302  };
303 
304  template<>
305  struct NormTraits<FFTWComplex>
306  {
307  typedef FFTWComplex Type;
308  typedef fftw_real SquaredNormType;
309  typedef fftw_real NormType;
310  };
311 
312  template <>
313  struct PromoteTraits<fftw_complex, fftw_complex>
314  {
315  typedef fftw_complex Promote;
316  };
317 
318  template <>
319  struct PromoteTraits<fftw_complex, double>
320  {
321  typedef fftw_complex Promote;
322  };
323 
324  template <>
325  struct PromoteTraits<double, fftw_complex>
326  {
327  typedef fftw_complex Promote;
328  };
329 
330  template <>
331  struct PromoteTraits<FFTWComplex, FFTWComplex>
332  {
333  typedef FFTWComplex Promote;
334  };
335 
336  template <>
337  struct PromoteTraits<FFTWComplex, double>
338  {
339  typedef FFTWComplex Promote;
340  };
341 
342  template <>
343  struct PromoteTraits<double, FFTWComplex>
344  {
345  typedef FFTWComplex Promote;
346  };
347  \endcode
348 
349  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
350  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
351  Namespace: vigra
352 
353 */
354 template<>
355 struct NumericTraits<fftw_complex>
356 {
357  typedef fftw_complex Type;
358  typedef fftw_complex Promote;
359  typedef fftw_complex RealPromote;
360  typedef fftw_complex ComplexPromote;
361  typedef fftw_real ValueType;
362 
363  typedef VigraFalseType isIntegral;
364  typedef VigraFalseType isScalar;
365  typedef NumericTraits<fftw_real>::isSigned isSigned;
366  typedef VigraFalseType isOrdered;
367  typedef VigraTrueType isComplex;
368 
369  static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
370  static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
371  static FFTWComplex nonZero() { return one(); }
372 
373  static const Promote & toPromote(const Type & v) { return v; }
374  static const RealPromote & toRealPromote(const Type & v) { return v; }
375  static const Type & fromPromote(const Promote & v) { return v; }
376  static const Type & fromRealPromote(const RealPromote & v) { return v; }
377 };
378 
379 template<>
380 struct NumericTraits<FFTWComplex>
381 {
382  typedef FFTWComplex Type;
383  typedef FFTWComplex Promote;
384  typedef FFTWComplex RealPromote;
385  typedef FFTWComplex ComplexPromote;
386  typedef fftw_real ValueType;
387 
388  typedef VigraFalseType isIntegral;
389  typedef VigraFalseType isScalar;
390  typedef NumericTraits<fftw_real>::isSigned isSigned;
391  typedef VigraFalseType isOrdered;
392  typedef VigraTrueType isComplex;
393 
394  static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
395  static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
396  static FFTWComplex nonZero() { return one(); }
397 
398  static const Promote & toPromote(const Type & v) { return v; }
399  static const RealPromote & toRealPromote(const Type & v) { return v; }
400  static const Type & fromPromote(const Promote & v) { return v; }
401  static const Type & fromRealPromote(const RealPromote & v) { return v; }
402 };
403 
404 template<>
405 struct NormTraits<fftw_complex>
406 {
407  typedef fftw_complex Type;
408  typedef fftw_real SquaredNormType;
409  typedef fftw_real NormType;
410 };
411 
412 template<>
413 struct NormTraits<FFTWComplex>
414 {
415  typedef FFTWComplex Type;
416  typedef fftw_real SquaredNormType;
417  typedef fftw_real NormType;
418 };
419 
420 template <>
421 struct PromoteTraits<fftw_complex, fftw_complex>
422 {
423  typedef fftw_complex Promote;
424 };
425 
426 template <>
427 struct PromoteTraits<fftw_complex, double>
428 {
429  typedef fftw_complex Promote;
430 };
431 
432 template <>
433 struct PromoteTraits<double, fftw_complex>
434 {
435  typedef fftw_complex Promote;
436 };
437 
438 template <>
439 struct PromoteTraits<FFTWComplex, FFTWComplex>
440 {
441  typedef FFTWComplex Promote;
442 };
443 
444 template <>
445 struct PromoteTraits<FFTWComplex, double>
446 {
447  typedef FFTWComplex Promote;
448 };
449 
450 template <>
451 struct PromoteTraits<double, FFTWComplex>
452 {
453  typedef FFTWComplex Promote;
454 };
455 
456 template<class T>
457 struct CanSkipInitialization<std::complex<T> >
458 {
459  typedef typename CanSkipInitialization<T>::type type;
460  static const bool value = type::asBool;
461 };
462 
463 template<>
464 struct CanSkipInitialization<FFTWComplex>
465 {
466  typedef CanSkipInitialization<fftw_real>::type type;
467  static const bool value = type::asBool;
468 };
469 
470 
471 
472 /********************************************************/
473 /* */
474 /* FFTWComplex Operations */
475 /* */
476 /********************************************************/
477 
478 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex
479 
480  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
481  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
482 
483  These functions fulfill the requirements of an Algebraic Field.
484  Return types are determined according to \ref FFTWComplexTraits.
485 
486  Namespace: vigra
487  <p>
488 
489  */
490 //@{
491  /// equal
492 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) {
493  return a.re() == b.re() && a.im() == b.im();
494 }
495 
496  /// not equal
497 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) {
498  return a.re() != b.re() || a.im() != b.im();
499 }
500 
501  /// add-assignment
503  a.re() += b.re();
504  a.im() += b.im();
505  return a;
506 }
507 
508  /// subtract-assignment
510  a.re() -= b.re();
511  a.im() -= b.im();
512  return a;
513 }
514 
515  /// multiply-assignment
517  FFTWComplex::value_type t = a.re()*b.re()-a.im()*b.im();
518  a.im() = a.re()*b.im()+a.im()*b.re();
519  a.re() = t;
520  return a;
521 }
522 
523  /// divide-assignment
526  FFTWComplex::value_type t = (a.re()*b.re()+a.im()*b.im())/sm;
527  a.im() = (b.re()*a.im()-a.re()*b.im())/sm;
528  a.re() = t;
529  return a;
530 }
531 
532  /// multiply-assignment with scalar double
533 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) {
534  a.re() *= b;
535  a.im() *= b;
536  return a;
537 }
538 
539  /// divide-assignment with scalar double
540 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) {
541  a.re() /= b;
542  a.im() /= b;
543  return a;
544 }
545 
546  /// addition
548  a += b;
549  return a;
550 }
551 
552  /// subtraction
554  a -= b;
555  return a;
556 }
557 
558  /// multiplication
559 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) {
560  a *= b;
561  return a;
562 }
563 
564  /// right multiplication with scalar double
565 inline FFTWComplex operator *(FFTWComplex a, const double &b) {
566  a *= b;
567  return a;
568 }
569 
570  /// left multiplication with scalar double
571 inline FFTWComplex operator *(const double &a, FFTWComplex b) {
572  b *= a;
573  return b;
574 }
575 
576  /// division
577 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) {
578  a /= b;
579  return a;
580 }
581 
582  /// right division with scalar double
583 inline FFTWComplex operator /(FFTWComplex a, const double &b) {
584  a /= b;
585  return a;
586 }
587 
588 using VIGRA_CSTD::abs;
589 
590  /// absolute value (= magnitude)
592 {
593  return a.magnitude();
594 }
595 
596  /// complex conjugate
597 inline FFTWComplex conj(const FFTWComplex &a)
598 {
599  return FFTWComplex(a.re(), -a.im());
600 }
601 
602  /// norm (= magnitude)
604 {
605  return a.magnitude();
606 }
607 
608  /// squared norm (= squared magnitude)
610 {
611  return a.squaredMagnitude();
612 }
613 
614 //@}
615 
616 
617 /** \addtogroup StandardImageTypes
618 */
619 //@{
620 
621 /********************************************************/
622 /* */
623 /* FFTWRealImage */
624 /* */
625 /********************************************************/
626 
627  /** Float (<tt>fftw_real</tt>) image.
628 
629  The type <tt>fftw_real</tt> is defined as <tt>double</tt> (in FFTW 2 it used to be
630  either <tt>float</tt> or <tt>double</tt>, as specified during compilation of FFTW).
631  FFTWRealImage uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
632  their const counterparts to access the data.
633 
634  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
635  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
636  Namespace: vigra
637  */
639 
640 /********************************************************/
641 /* */
642 /* FFTWComplexImage */
643 /* */
644 /********************************************************/
645 
646 template<>
647 struct IteratorTraits<
649 {
651  typedef Iterator iterator;
652  typedef BasicImageIterator<FFTWComplex, FFTWComplex **> mutable_iterator;
654  typedef iterator::iterator_category iterator_category;
655  typedef iterator::value_type value_type;
656  typedef iterator::reference reference;
657  typedef iterator::index_reference index_reference;
658  typedef iterator::pointer pointer;
659  typedef iterator::difference_type difference_type;
660  typedef iterator::row_iterator row_iterator;
661  typedef iterator::column_iterator column_iterator;
662  typedef VectorAccessor<FFTWComplex> default_accessor;
663  typedef VectorAccessor<FFTWComplex> DefaultAccessor;
664  typedef VigraTrueType hasConstantStrides;
665 };
666 
667 template<>
668 struct IteratorTraits<
669  ConstBasicImageIterator<FFTWComplex, FFTWComplex **> >
670 {
671  typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> Iterator;
672  typedef Iterator iterator;
673  typedef BasicImageIterator<FFTWComplex, FFTWComplex **> mutable_iterator;
674  typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> const_iterator;
675  typedef iterator::iterator_category iterator_category;
676  typedef iterator::value_type value_type;
677  typedef iterator::reference reference;
678  typedef iterator::index_reference index_reference;
679  typedef iterator::pointer pointer;
680  typedef iterator::difference_type difference_type;
681  typedef iterator::row_iterator row_iterator;
682  typedef iterator::column_iterator column_iterator;
683  typedef VectorAccessor<FFTWComplex> default_accessor;
684  typedef VectorAccessor<FFTWComplex> DefaultAccessor;
685  typedef VigraTrueType hasConstantStrides;
686 };
687 
688  /** Complex (FFTWComplex) image.
689  It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
690  their const counterparts to access the data.
691 
692  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
693  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
694  Namespace: vigra
695  */
697 
698 //@}
699 
700 /********************************************************/
701 /* */
702 /* FFTWComplex-Accessors */
703 /* */
704 /********************************************************/
705 
706 /** \addtogroup DataAccessors
707 */
708 //@{
709 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex
710 
711  Encapsulate access to pixels of type FFTWComplex
712 */
713 //@{
714  /** Encapsulate access to the the real part of a complex number.
715 
716  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
717  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
718  Namespace: vigra
719  */
721 {
722  public:
723 
724  /// The accessor's value type.
725  typedef fftw_real value_type;
726 
727  /// Read real part at iterator position.
728  template <class ITERATOR>
729  value_type operator()(ITERATOR const & i) const {
730  return (*i).re();
731  }
732 
733  /// Read real part at offset from iterator position.
734  template <class ITERATOR, class DIFFERENCE>
735  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
736  return i[d].re();
737  }
738 
739  /// Write real part at iterator position from a scalar.
740  template <class ITERATOR>
741  void set(value_type const & v, ITERATOR const & i) const {
742  (*i).re()= v;
743  }
744 
745  /// Write real part at offset from iterator position from a scalar.
746  template <class ITERATOR, class DIFFERENCE>
747  void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
748  i[d].re()= v;
749  }
750 
751  /// Write real part at iterator position into a scalar.
752  template <class ITERATOR>
753  void set(FFTWComplex const & v, ITERATOR const & i) const {
754  *i = v.re();
755  }
756 
757  /// Write real part at offset from iterator position into a scalar.
758  template <class ITERATOR, class DIFFERENCE>
759  void set(FFTWComplex const & v, ITERATOR const & i, DIFFERENCE d) const {
760  i[d] = v.re();
761  }
762 };
763 
764  /** Encapsulate access to the the imaginary part of a complex number.
765 
766  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
767  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
768  Namespace: vigra
769  */
771 {
772  public:
773  /// The accessor's value type.
774  typedef fftw_real value_type;
775 
776  /// Read imaginary part at iterator position.
777  template <class ITERATOR>
778  value_type operator()(ITERATOR const & i) const {
779  return (*i).im();
780  }
781 
782  /// Read imaginary part at offset from iterator position.
783  template <class ITERATOR, class DIFFERENCE>
784  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
785  return i[d].im();
786  }
787 
788  /// Write imaginary part at iterator position from a scalar.
789  template <class ITERATOR>
790  void set(value_type const & v, ITERATOR const & i) const {
791  (*i).im()= v;
792  }
793 
794  /// Write imaginary part at offset from iterator position from a scalar.
795  template <class ITERATOR, class DIFFERENCE>
796  void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
797  i[d].im()= v;
798  }
799 
800  /// Write imaginary part at iterator position into a scalar.
801  template <class ITERATOR>
802  void set(FFTWComplex const & v, ITERATOR const & i) const {
803  *i = v.im();
804  }
805 
806  /// Write imaginary part at offset from iterator position into a scalar.
807  template <class ITERATOR, class DIFFERENCE>
808  void set(FFTWComplex const & v, ITERATOR const & i, DIFFERENCE d) const {
809  i[d] = v.im();
810  }
811 };
812 
813  /** Write a real number into a complex one. The imaginary part is set
814  to 0.
815 
816  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
817  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
818  Namespace: vigra
819  */
821 {
822  public:
823  /// The accessor's value type.
824  typedef fftw_real value_type;
825 
826  /** Write real number at iterator position. Set imaginary part
827  to 0.
828  */
829  template <class ITERATOR>
830  void set(value_type const & v, ITERATOR const & i) const {
831  (*i).re()= v;
832  (*i).im()= 0;
833  }
834 
835  /** Write real number at offset from iterator position. Set imaginary part
836  to 0.
837  */
838  template <class ITERATOR, class DIFFERENCE>
839  void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
840  i[d].re()= v;
841  i[d].im()= 0;
842  }
843 };
844 
845  /** Calculate magnitude of complex number on the fly.
846 
847  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
848  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
849  Namespace: vigra
850  */
852 {
853  public:
854  /// The accessor's value type.
855  typedef fftw_real value_type;
856 
857  /// Read magnitude at iterator position.
858  template <class ITERATOR>
859  value_type operator()(ITERATOR const & i) const {
860  return (*i).magnitude();
861  }
862 
863  /// Read magnitude at offset from iterator position.
864  template <class ITERATOR, class DIFFERENCE>
865  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
866  return (i[d]).magnitude();
867  }
868 };
869 
870  /** Calculate phase of complex number on the fly.
871 
872  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br>
873  <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br>
874  Namespace: vigra
875  */
877 {
878  public:
879  /// The accessor's value type.
880  typedef fftw_real value_type;
881 
882  /// Read phase at iterator position.
883  template <class ITERATOR>
884  value_type operator()(ITERATOR const & i) const {
885  return (*i).phase();
886  }
887 
888  /// Read phase at offset from iterator position.
889  template <class ITERATOR, class DIFFERENCE>
890  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
891  return (i[d]).phase();
892  }
893 };
894 
895 //@}
896 //@}
897 
898 /********************************************************/
899 /* */
900 /* Fourier Transform */
901 /* */
902 /********************************************************/
903 
904 /** \addtogroup FourierTransform Fast Fourier Transform
905 
906  This documentation describes the VIGRA interface to FFTW version 3. The interface
907  to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated.
908 
909  VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
910  Transform</a> package to perform Fourier transformations. VIGRA
911  provides a wrapper for FFTW's complex number type (FFTWComplex),
912  but FFTW's functions are used verbatim. If the image is stored as
913  a FFTWComplexImage, the simplest call to an FFT function is like this:
914 
915  \code
916  vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
917  ... // fill image with data
918 
919  // create a plan with estimated performance optimization
920  fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
921  (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
922  FFTW_FORWARD, FFTW_ESTIMATE );
923  // calculate FFT (this can be repeated as often as needed,
924  // with fresh data written into the source array)
925  fftw_execute(forwardPlan);
926 
927  // release the plan memory
928  fftw_destroy_plan(forwardPlan);
929 
930  // likewise for the inverse transform
931  fftw_plan backwardPlan = fftw_plan_dft_2d(height, width,
932  (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(),
933  FFTW_BACKWARD, FFTW_ESTIMATE);
934  fftw_execute(backwardPlan);
935  fftw_destroy_plan(backwardPlan);
936 
937  // do not forget to normalize the result according to the image size
938  transformImage(srcImageRange(spatial), destImage(spatial),
939  std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height));
940  \endcode
941 
942  Note that in the creation of a plan, the height must be given
943  first. Note also that <TT>spatial.begin()</TT> may only be passed
944  to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the
945  entire image. When you want to restrict operation to an ROI, you
946  can create a copy of the ROI in an image of appropriate size, or
947  you may use the Guru interface to FFTW.
948 
949  More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
950 
951  FFTW produces fourier images that have the DC component (the
952  origin of the Fourier space) in the upper left corner. Often, one
953  wants the origin in the center of the image, so that frequencies
954  always increase towards the border of the image. This can be
955  achieved by calling \ref moveDCToCenter(). The inverse
956  transformation is done by \ref moveDCToUpperLeft().
957 
958  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br>
959  Namespace: vigra
960 */
961 
962 /** \addtogroup FourierTransform
963 */
964 //@{
965 
966 /********************************************************/
967 /* */
968 /* moveDCToCenter */
969 /* */
970 /********************************************************/
971 
972 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
973  in the image center.
974 
975  FFTW produces fourier images where the DC component (origin of
976  fourier space) is located in the upper left corner of the
977  image. The quadrants are placed like this (using a 4x4 image for
978  example):
979 
980  \code
981  DC 4 3 3
982  4 4 3 3
983  1 1 2 2
984  1 1 2 2
985  \endcode
986 
987  After applying the function, the quadrants are at their usual places:
988 
989  \code
990  2 2 1 1
991  2 2 1 1
992  3 3 DC 4
993  3 3 4 4
994  \endcode
995 
996  This transformation can be reversed by \ref moveDCToUpperLeft().
997  Note that the transformation must not be executed in place - input
998  and output images must be different.
999 
1000  <b> Declarations:</b>
1001 
1002  pass arguments explicitly:
1003  \code
1004  namespace vigra {
1005  template <class SrcImageIterator, class SrcAccessor,
1006  class DestImageIterator, class DestAccessor>
1007  void moveDCToCenter(SrcImageIterator src_upperleft,
1008  SrcImageIterator src_lowerright, SrcAccessor sa,
1009  DestImageIterator dest_upperleft, DestAccessor da);
1010  }
1011  \endcode
1012 
1013 
1014  use argument objects in conjunction with \ref ArgumentObjectFactories :
1015  \code
1016  namespace vigra {
1017  template <class SrcImageIterator, class SrcAccessor,
1018  class DestImageIterator, class DestAccessor>
1019  void moveDCToCenter(
1020  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1021  pair<DestImageIterator, DestAccessor> dest);
1022  }
1023  \endcode
1024 
1025  <b> Usage:</b>
1026 
1027  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br>
1028  Namespace: vigra
1029 
1030  \code
1031  vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
1032  ... // fill image with data
1033 
1034  // create a plan with estimated performance optimization
1035  fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
1036  (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
1037  FFTW_FORWARD, FFTW_ESTIMATE );
1038  // calculate FFT
1039  fftw_execute(forwardPlan);
1040 
1041  vigra::FFTWComplexImage rearrangedFourier(width, height);
1042  moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
1043 
1044  // delete the plan
1045  fftw_destroy_plan(forwardPlan);
1046  \endcode
1047 */
1048 doxygen_overloaded_function(template <...> void moveDCToCenter)
1049 
1050 template <class SrcImageIterator, class SrcAccessor,
1051  class DestImageIterator, class DestAccessor>
1052 void moveDCToCenter(SrcImageIterator src_upperleft,
1053  SrcImageIterator src_lowerright, SrcAccessor sa,
1054  DestImageIterator dest_upperleft, DestAccessor da)
1055 {
1056  int w = int(src_lowerright.x - src_upperleft.x);
1057  int h = int(src_lowerright.y - src_upperleft.y);
1058  int w1 = w/2;
1059  int h1 = h/2;
1060  int w2 = (w+1)/2;
1061  int h2 = (h+1)/2;
1062 
1063  // 2. Quadrant zum 4.
1064  copyImage(srcIterRange(src_upperleft,
1065  src_upperleft + Diff2D(w2, h2), sa),
1066  destIter (dest_upperleft + Diff2D(w1, h1), da));
1067 
1068  // 4. Quadrant zum 2.
1069  copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1070  src_lowerright, sa),
1071  destIter (dest_upperleft, da));
1072 
1073  // 1. Quadrant zum 3.
1074  copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1075  src_upperleft + Diff2D(w, h2), sa),
1076  destIter (dest_upperleft + Diff2D(0, h1), da));
1077 
1078  // 3. Quadrant zum 1.
1079  copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1080  src_upperleft + Diff2D(w2, h), sa),
1081  destIter (dest_upperleft + Diff2D(w1, 0), da));
1082 }
1083 
1084 template <class SrcImageIterator, class SrcAccessor,
1085  class DestImageIterator, class DestAccessor>
1086 inline void moveDCToCenter(
1087  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1088  pair<DestImageIterator, DestAccessor> dest)
1089 {
1090  moveDCToCenter(src.first, src.second, src.third,
1091  dest.first, dest.second);
1092 }
1093 
1094 /********************************************************/
1095 /* */
1096 /* moveDCToUpperLeft */
1097 /* */
1098 /********************************************************/
1099 
1100 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
1101  in the image's upper left.
1102 
1103  This function is the inversion of \ref moveDCToCenter(). See there
1104  for declarations and a usage example.
1105 
1106  <b> Declarations:</b>
1107 
1108  pass arguments explicitly:
1109  \code
1110  namespace vigra {
1111  template <class SrcImageIterator, class SrcAccessor,
1112  class DestImageIterator, class DestAccessor>
1113  void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1114  SrcImageIterator src_lowerright, SrcAccessor sa,
1115  DestImageIterator dest_upperleft, DestAccessor da);
1116  }
1117  \endcode
1118 
1119 
1120  use argument objects in conjunction with \ref ArgumentObjectFactories :
1121  \code
1122  namespace vigra {
1123  template <class SrcImageIterator, class SrcAccessor,
1124  class DestImageIterator, class DestAccessor>
1125  void moveDCToUpperLeft(
1126  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1127  pair<DestImageIterator, DestAccessor> dest);
1128  }
1129  \endcode
1130 */
1132 
1133 template <class SrcImageIterator, class SrcAccessor,
1134  class DestImageIterator, class DestAccessor>
1135 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1136  SrcImageIterator src_lowerright, SrcAccessor sa,
1137  DestImageIterator dest_upperleft, DestAccessor da)
1138 {
1139  int w = int(src_lowerright.x - src_upperleft.x);
1140  int h = int(src_lowerright.y - src_upperleft.y);
1141  int w2 = w/2;
1142  int h2 = h/2;
1143  int w1 = (w+1)/2;
1144  int h1 = (h+1)/2;
1145 
1146  // 2. Quadrant zum 4.
1147  copyImage(srcIterRange(src_upperleft,
1148  src_upperleft + Diff2D(w2, h2), sa),
1149  destIter (dest_upperleft + Diff2D(w1, h1), da));
1150 
1151  // 4. Quadrant zum 2.
1152  copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1153  src_lowerright, sa),
1154  destIter (dest_upperleft, da));
1155 
1156  // 1. Quadrant zum 3.
1157  copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1158  src_upperleft + Diff2D(w, h2), sa),
1159  destIter (dest_upperleft + Diff2D(0, h1), da));
1160 
1161  // 3. Quadrant zum 1.
1162  copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1163  src_upperleft + Diff2D(w2, h), sa),
1164  destIter (dest_upperleft + Diff2D(w1, 0), da));
1165 }
1166 
1167 template <class SrcImageIterator, class SrcAccessor,
1168  class DestImageIterator, class DestAccessor>
1169 inline void moveDCToUpperLeft(
1170  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1171  pair<DestImageIterator, DestAccessor> dest)
1172 {
1173  moveDCToUpperLeft(src.first, src.second, src.third,
1174  dest.first, dest.second);
1175 }
1176 
1177 template <class DestImageIterator, class DestAccessor>
1178 void fftShift(DestImageIterator upperleft,
1179  DestImageIterator lowerright, DestAccessor da)
1180 {
1181  int w = int(lowerright.x - upperleft.x);
1182  int h = int(lowerright.y - upperleft.y);
1183  int w2 = w/2;
1184  int h2 = h/2;
1185  int w1 = (w+1)/2;
1186  int h1 = (h+1)/2;
1187 
1188  // 2. Quadrant zum 4.
1189  swapImageData(destIterRange(upperleft,
1190  upperleft + Diff2D(w2, h2), da),
1191  destIter (upperleft + Diff2D(w1, h1), da));
1192 
1193  // 1. Quadrant zum 3.
1194  swapImageData(destIterRange(upperleft + Diff2D(w2, 0),
1195  upperleft + Diff2D(w, h2), da),
1196  destIter (upperleft + Diff2D(0, h1), da));
1197 }
1198 
1199 template <class DestImageIterator, class DestAccessor>
1200 inline void fftShift(
1201  triple<DestImageIterator, DestImageIterator, DestAccessor> dest)
1202 {
1203  fftShift(dest.first, dest.second, dest.third);
1204 }
1205 
1206 
1207 namespace detail {
1208 
1209 template <class T>
1210 void
1211 fourierTransformImpl(FFTWComplexImage::const_traverser sul,
1214 {
1215  int w = int(slr.x - sul.x);
1216  int h = int(slr.y - sul.y);
1217 
1218  FFTWComplexImage sworkImage, dworkImage;
1219 
1220  fftw_complex * srcPtr = (fftw_complex *)(&*sul);
1221  fftw_complex * destPtr = (fftw_complex *)(&*dul);
1222 
1223  // test for right memory layout (fftw expects a 2*width*height floats array)
1224  if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1))))
1225  {
1226  sworkImage.resize(w, h);
1227  copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
1228  srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft()));
1229  }
1230  if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1231  {
1232  dworkImage.resize(w, h);
1233  destPtr = (fftw_complex *)(&(*dworkImage.upperLeft()));
1234  }
1235 
1236  fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
1237  fftw_execute(plan);
1238  fftw_destroy_plan(plan);
1239 
1240  if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1241  {
1242  copyImage(srcImageRange(dworkImage), destIter(dul, dest));
1243  }
1244 }
1245 
1246 } // namespace detail
1247 
1248 /********************************************************/
1249 /* */
1250 /* fourierTransform */
1251 /* */
1252 /********************************************************/
1253 
1254 /** \brief Compute forward and inverse Fourier transforms.
1255 
1256  In the forward direction, the input image may be scalar or complex, and the output image
1257  is always complex. In the inverse direction, both input and output must be complex.
1258 
1259  <b> Declarations:</b>
1260 
1261  pass arguments explicitly:
1262  \code
1263  namespace vigra {
1264  template <class SrcImageIterator, class SrcAccessor>
1265  void fourierTransform(SrcImageIterator srcUpperLeft,
1266  SrcImageIterator srcLowerRight, SrcAccessor src,
1267  FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest);
1268 
1269  void
1270  fourierTransformInverse(FFTWComplexImage::const_traverser sul,
1271  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1272  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
1273  }
1274  \endcode
1275 
1276  use argument objects in conjunction with \ref ArgumentObjectFactories :
1277  \code
1278  namespace vigra {
1279  template <class SrcImageIterator, class SrcAccessor>
1280  void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1281  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
1282 
1283  void
1284  fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
1285  FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
1286  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
1287  }
1288  \endcode
1289 
1290  <b> Usage:</b>
1291 
1292  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br>
1293  Namespace: vigra
1294 
1295  \code
1296  // compute complex Fourier transform of a real image
1297  vigra::DImage src(w, h);
1298  vigra::FFTWComplexImage fourier(w, h);
1299 
1300  fourierTransform(srcImageRange(src), destImage(fourier));
1301 
1302  // compute inverse Fourier transform
1303  // note that both source and destination image must be of type vigra::FFTWComplexImage
1304  vigra::FFTWComplexImage inverseFourier(w, h);
1305 
1306  fourierTransform(srcImageRange(fourier), destImage(inverseFourier));
1307  \endcode
1308 */
1309 doxygen_overloaded_function(template <...> void fourierTransform)
1310 
1311 inline void
1315 {
1316  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
1317 }
1318 
1319 template <class SrcImageIterator, class SrcAccessor>
1320 void fourierTransform(SrcImageIterator srcUpperLeft,
1321  SrcImageIterator srcLowerRight, SrcAccessor sa,
1323 {
1324  // copy real input images into a complex one...
1325  int w= srcLowerRight.x - srcUpperLeft.x;
1326  int h= srcLowerRight.y - srcUpperLeft.y;
1327 
1328  FFTWComplexImage workImage(w, h);
1329  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1330  destImage(workImage, FFTWWriteRealAccessor()));
1331 
1332  // ...and call the complex -> complex version of the algorithm
1333  FFTWComplexImage const & cworkImage = workImage;
1334  fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1335  destUpperLeft, da);
1336 }
1337 
1338 template <class SrcImageIterator, class SrcAccessor>
1339 inline
1340 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1341  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
1342 {
1343  fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
1344 }
1345 
1346 /** \brief Compute inverse Fourier transforms.
1347 
1348  See \ref fourierTransform() for details.
1349 */
1350 inline void
1354 {
1355  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
1356 }
1357 
1358 template <class DestImageIterator, class DestAccessor>
1361  DestImageIterator dul, DestAccessor dest)
1362 {
1363  int w = slr.x - sul.x;
1364  int h = slr.y - sul.y;
1365 
1366  FFTWComplexImage workImage(w, h);
1367  fourierTransformInverse(sul, slr, src, workImage.upperLeft(), workImage.accessor());
1368  copyImage(srcImageRange(workImage), destIter(dul, dest));
1369 }
1370 
1371 
1372 template <class DestImageIterator, class DestAccessor>
1373 inline void
1376  pair<DestImageIterator, DestAccessor> dest)
1377 {
1378  fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second);
1379 }
1380 
1381 
1382 
1383 /********************************************************/
1384 /* */
1385 /* applyFourierFilter */
1386 /* */
1387 /********************************************************/
1388 
1389 /** \brief Apply a filter (defined in the frequency domain) to an image.
1390 
1391  After transferring the image into the frequency domain, it is
1392  multiplied pixel-wise with the filter and transformed back. The
1393  result is put into the given destination image which must have the right size.
1394  The result will be normalized to compensate for the two FFTs.
1395 
1396  If the destination image is scalar, only the real part of the result image is
1397  retained. In this case, you are responsible for choosing a filter image
1398  which ensures a zero imaginary part of the result (e.g. use a real, even symmetric
1399  filter image, or a purely imaginary, odd symmetric on).
1400 
1401  The DC entry of the filter must be in the upper left, which is the
1402  position where FFTW expects it (see \ref moveDCToUpperLeft()).
1403 
1404  <b> Declarations:</b>
1405 
1406  pass arguments explicitly:
1407  \code
1408  namespace vigra {
1409  template <class SrcImageIterator, class SrcAccessor,
1410  class FilterImageIterator, class FilterAccessor,
1411  class DestImageIterator, class DestAccessor>
1412  void applyFourierFilter(SrcImageIterator srcUpperLeft,
1413  SrcImageIterator srcLowerRight, SrcAccessor sa,
1414  FilterImageIterator filterUpperLeft, FilterAccessor fa,
1415  DestImageIterator destUpperLeft, DestAccessor da);
1416  }
1417  \endcode
1418 
1419  use argument objects in conjunction with \ref ArgumentObjectFactories :
1420  \code
1421  namespace vigra {
1422  template <class SrcImageIterator, class SrcAccessor,
1423  class FilterImageIterator, class FilterAccessor,
1424  class DestImageIterator, class DestAccessor>
1425  void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1426  pair<FilterImageIterator, FilterAccessor> filter,
1427  pair<DestImageIterator, DestAccessor> dest);
1428  }
1429  \endcode
1430 
1431  <b> Usage:</b>
1432 
1433  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br>
1434  Namespace: vigra
1435 
1436  \code
1437  // create a Gaussian filter in Fourier space
1438  vigra::FImage gaussFilter(w, h), filter(w, h);
1439  for(int y=0; y<h; ++y)
1440  for(int x=0; x<w; ++x)
1441  {
1442  xx = float(x - w / 2) / w;
1443  yy = float(y - h / 2) / h;
1444 
1445  gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
1446  }
1447 
1448  // applyFourierFilter() expects the filter's DC in the upper left
1449  moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
1450 
1451  vigra::FFTWComplexImage result(w, h);
1452 
1453  vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
1454  \endcode
1455 
1456  For inspection of the result, \ref FFTWMagnitudeAccessor might be
1457  useful. If you want to apply the same filter repeatedly, it may be more
1458  efficient to use the FFTW functions directly with FFTW plans optimized
1459  for good performance.
1460 */
1462 
1463 template <class SrcImageIterator, class SrcAccessor,
1464  class FilterImageIterator, class FilterAccessor,
1465  class DestImageIterator, class DestAccessor>
1466 void applyFourierFilter(SrcImageIterator srcUpperLeft,
1467  SrcImageIterator srcLowerRight, SrcAccessor sa,
1468  FilterImageIterator filterUpperLeft, FilterAccessor fa,
1469  DestImageIterator destUpperLeft, DestAccessor da)
1470 {
1471  // copy real input images into a complex one...
1472  int w = int(srcLowerRight.x - srcUpperLeft.x);
1473  int h = int(srcLowerRight.y - srcUpperLeft.y);
1474 
1475  FFTWComplexImage workImage(w, h);
1476  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1477  destImage(workImage, FFTWWriteRealAccessor()));
1478 
1479  // ...and call the impl
1480  FFTWComplexImage const & cworkImage = workImage;
1481  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1482  filterUpperLeft, fa,
1483  destUpperLeft, da);
1484 }
1485 
1486 template <class FilterImageIterator, class FilterAccessor,
1487  class DestImageIterator, class DestAccessor>
1488 inline
1489 void applyFourierFilter(
1490  FFTWComplexImage::const_traverser srcUpperLeft,
1491  FFTWComplexImage::const_traverser srcLowerRight,
1493  FilterImageIterator filterUpperLeft, FilterAccessor fa,
1494  DestImageIterator destUpperLeft, DestAccessor da)
1495 {
1496  int w = srcLowerRight.x - srcUpperLeft.x;
1497  int h = srcLowerRight.y - srcUpperLeft.y;
1498 
1499  // test for right memory layout (fftw expects a 2*width*height floats array)
1500  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
1501  applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
1502  filterUpperLeft, fa,
1503  destUpperLeft, da);
1504  else
1505  {
1506  FFTWComplexImage workImage(w, h);
1507  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1508  destImage(workImage));
1509 
1510  FFTWComplexImage const & cworkImage = workImage;
1511  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1512  filterUpperLeft, fa,
1513  destUpperLeft, da);
1514  }
1515 }
1516 
1517 template <class SrcImageIterator, class SrcAccessor,
1518  class FilterImageIterator, class FilterAccessor,
1519  class DestImageIterator, class DestAccessor>
1520 inline
1521 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1522  pair<FilterImageIterator, FilterAccessor> filter,
1523  pair<DestImageIterator, DestAccessor> dest)
1524 {
1525  applyFourierFilter(src.first, src.second, src.third,
1526  filter.first, filter.second,
1527  dest.first, dest.second);
1528 }
1529 
1530 template <class FilterImageIterator, class FilterAccessor,
1531  class DestImageIterator, class DestAccessor>
1532 void applyFourierFilterImpl(
1533  FFTWComplexImage::const_traverser srcUpperLeft,
1534  FFTWComplexImage::const_traverser srcLowerRight,
1536  FilterImageIterator filterUpperLeft, FilterAccessor fa,
1537  DestImageIterator destUpperLeft, DestAccessor da)
1538 {
1539  int w = int(srcLowerRight.x - srcUpperLeft.x);
1540  int h = int(srcLowerRight.y - srcUpperLeft.y);
1541 
1542  FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
1543 
1544  // FFT from srcImage to complexResultImg
1545  fftw_plan forwardPlan=
1546  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
1547  (fftw_complex *)complexResultImg.begin(),
1548  FFTW_FORWARD, FFTW_ESTIMATE );
1549  fftw_execute(forwardPlan);
1550  fftw_destroy_plan(forwardPlan);
1551 
1552  // convolve in freq. domain (in complexResultImg)
1553  combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
1554  destImage(complexResultImg), std::multiplies<FFTWComplex>());
1555 
1556  // FFT back into spatial domain (inplace in complexResultImg)
1557  fftw_plan backwardPlan=
1558  fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
1559  (fftw_complex *)complexResultImg.begin(),
1560  FFTW_BACKWARD, FFTW_ESTIMATE);
1561  fftw_execute(backwardPlan);
1562  fftw_destroy_plan(backwardPlan);
1563 
1564  typedef typename
1565  NumericTraits<typename DestAccessor::value_type>::isScalar
1566  isScalarResult;
1567 
1568  // normalization (after FFTs), maybe stripping imaginary part
1569  applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
1570  isScalarResult());
1571 }
1572 
1573 template <class DestImageIterator, class DestAccessor>
1574 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
1575  DestImageIterator destUpperLeft,
1576  DestAccessor da,
1577  VigraFalseType)
1578 {
1579  double normFactor= 1.0/(srcImage.width() * srcImage.height());
1580 
1581  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
1582  {
1583  DestImageIterator dIt= destUpperLeft;
1584  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
1585  {
1586  da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
1587  da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
1588  }
1589  }
1590 }
1591 
1592 inline
1593 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
1594  FFTWComplexImage::traverser destUpperLeft,
1596  VigraFalseType)
1597 {
1598  transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
1599  linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height())));
1600 }
1601 
1602 template <class DestImageIterator, class DestAccessor>
1603 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
1604  DestImageIterator destUpperLeft,
1605  DestAccessor da,
1606  VigraTrueType)
1607 {
1608  double normFactor= 1.0/(srcImage.width() * srcImage.height());
1609 
1610  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
1611  {
1612  DestImageIterator dIt= destUpperLeft;
1613  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
1614  da.set(srcImage(x, y).re()*normFactor, dIt);
1615  }
1616 }
1617 
1618 /**********************************************************/
1619 /* */
1620 /* applyFourierFilterFamily */
1621 /* */
1622 /**********************************************************/
1623 
1624 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
1625 
1626  This provides the same functionality as \ref applyFourierFilter(),
1627  but applying several filters at once allows to avoid
1628  repeated Fourier transforms of the source image.
1629 
1630  Filters and result images must be stored in \ref vigra::ImageArray data
1631  structures. In contrast to \ref applyFourierFilter(), this function adjusts
1632  the size of the result images and the the length of the array.
1633 
1634  <b> Declarations:</b>
1635 
1636  pass arguments explicitly:
1637  \code
1638  namespace vigra {
1639  template <class SrcImageIterator, class SrcAccessor, class FilterType>
1640  void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
1641  SrcImageIterator srcLowerRight, SrcAccessor sa,
1642  const ImageArray<FilterType> &filters,
1643  ImageArray<FFTWComplexImage> &results)
1644  }
1645  \endcode
1646 
1647  use argument objects in conjunction with \ref ArgumentObjectFactories :
1648  \code
1649  namespace vigra {
1650  template <class SrcImageIterator, class SrcAccessor, class FilterType>
1651  void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1652  const ImageArray<FilterType> &filters,
1653  ImageArray<FFTWComplexImage> &results)
1654  }
1655  \endcode
1656 
1657  <b> Usage:</b>
1658 
1659  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br>
1660  Namespace: vigra
1661 
1662  \code
1663  // assuming the presence of a real-valued image named "image" to
1664  // be filtered in this example
1665 
1666  vigra::ImageArray<vigra::FImage> filters(16, image.size());
1667 
1668  for (int i=0; i<filters.size(); i++)
1669  // create some meaningful filters here
1670  createMyFilterOfScale(i, destImage(filters[i]));
1671 
1672  vigra::ImageArray<vigra::FFTWComplexImage> results();
1673 
1674  vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
1675  \endcode
1676 */
1678 
1679 template <class SrcImageIterator, class SrcAccessor,
1680  class FilterType, class DestImage>
1681 inline
1682 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1683  const ImageArray<FilterType> &filters,
1684  ImageArray<DestImage> &results)
1685 {
1686  applyFourierFilterFamily(src.first, src.second, src.third,
1687  filters, results);
1688 }
1689 
1690 template <class SrcImageIterator, class SrcAccessor,
1691  class FilterType, class DestImage>
1692 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
1693  SrcImageIterator srcLowerRight, SrcAccessor sa,
1694  const ImageArray<FilterType> &filters,
1695  ImageArray<DestImage> &results)
1696 {
1697  int w = int(srcLowerRight.x - srcUpperLeft.x);
1698  int h = int(srcLowerRight.y - srcUpperLeft.y);
1699 
1700  FFTWComplexImage workImage(w, h);
1701  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1702  destImage(workImage, FFTWWriteRealAccessor()));
1703 
1704  FFTWComplexImage const & cworkImage = workImage;
1705  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1706  filters, results);
1707 }
1708 
1709 template <class FilterType, class DestImage>
1710 inline
1712  FFTWComplexImage::const_traverser srcUpperLeft,
1713  FFTWComplexImage::const_traverser srcLowerRight,
1715  const ImageArray<FilterType> &filters,
1716  ImageArray<DestImage> &results)
1717 {
1718  int w= srcLowerRight.x - srcUpperLeft.x;
1719 
1720  // test for right memory layout (fftw expects a 2*width*height floats array)
1721  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
1722  applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
1723  filters, results);
1724  else
1725  {
1726  int h = srcLowerRight.y - srcUpperLeft.y;
1727  FFTWComplexImage workImage(w, h);
1728  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1729  destImage(workImage));
1730 
1731  FFTWComplexImage const & cworkImage = workImage;
1732  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1733  filters, results);
1734  }
1735 }
1736 
1737 template <class FilterType, class DestImage>
1738 void applyFourierFilterFamilyImpl(
1739  FFTWComplexImage::const_traverser srcUpperLeft,
1740  FFTWComplexImage::const_traverser srcLowerRight,
1742  const ImageArray<FilterType> &filters,
1743  ImageArray<DestImage> &results)
1744 {
1745  // FIXME: sa is not used
1746  // (maybe check if StandardAccessor, else copy?)
1747 
1748  // make sure the filter images have the right dimensions
1749  vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
1750  "applyFourierFilterFamily called with src image size != filters.imageSize()!");
1751 
1752  // make sure the result image array has the right dimensions
1753  results.resize(filters.size());
1754  results.resizeImages(filters.imageSize());
1755 
1756  // FFT from srcImage to freqImage
1757  int w = int(srcLowerRight.x - srcUpperLeft.x);
1758  int h = int(srcLowerRight.y - srcUpperLeft.y);
1759 
1760  FFTWComplexImage freqImage(w, h);
1761  FFTWComplexImage result(w, h);
1762 
1763  fftw_plan forwardPlan=
1764  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
1765  (fftw_complex *)freqImage.begin(),
1766  FFTW_FORWARD, FFTW_ESTIMATE );
1767  fftw_execute(forwardPlan);
1768  fftw_destroy_plan(forwardPlan);
1769 
1770  fftw_plan backwardPlan=
1771  fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
1772  (fftw_complex *)result.begin(),
1773  FFTW_BACKWARD, FFTW_ESTIMATE );
1774  typedef typename
1775  NumericTraits<typename DestImage::Accessor::value_type>::isScalar
1776  isScalarResult;
1777 
1778  // convolve with filters in freq. domain
1779  for (unsigned int i= 0; i < filters.size(); i++)
1780  {
1781  combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
1782  destImage(result), std::multiplies<FFTWComplex>());
1783 
1784  // FFT back into spatial domain (inplace in destImage)
1785  fftw_execute(backwardPlan);
1786 
1787  // normalization (after FFTs), maybe stripping imaginary part
1788  applyFourierFilterImplNormalization(result,
1789  results[i].upperLeft(), results[i].accessor(),
1790  isScalarResult());
1791  }
1792  fftw_destroy_plan(backwardPlan);
1793 }
1794 
1795 /********************************************************/
1796 /* */
1797 /* fourierTransformReal */
1798 /* */
1799 /********************************************************/
1800 
1801 /** \brief Real Fourier transforms for even and odd boundary conditions
1802  (aka. cosine and sine transforms).
1803 
1804 
1805  If the image is real and has even symmetry, its Fourier transform
1806  is also real and has even symmetry. The Fourier transform of a real image with odd
1807  symmetry is imaginary and has odd symmetry. In either case, only about a quarter
1808  of the pixels need to be stored because the rest can be calculated from the symmetry
1809  properties. This is especially useful, if the original image is implicitly assumed
1810  to have reflective or anti-reflective boundary conditions. Then the "negative"
1811  pixel locations are defined as
1812 
1813  \code
1814  even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1)
1815  odd (anti-reflective boundary conditions): f[-1] = 0
1816  f[-x] = -f[x-2] (x = 2,...,N-1)
1817  \endcode
1818 
1819  end similar at the other boundary (see the FFTW documentation for details).
1820  This has the advantage that more efficient Fourier transforms that use only
1821  real numbers can be implemented. These are also known as cosine and sine transforms
1822  respectively.
1823 
1824  If you use the odd transform it is important to note that in the Fourier domain,
1825  the DC component is always zero and is therefore dropped from the data structure.
1826  This means that index 0 in an odd symmetric Fourier domain image refers to
1827  the <i>first</i> harmonic. This is especially important if an image is first
1828  cosine transformed (even symmetry), then in the Fourier domain multiplied
1829  with an odd symmetric filter (e.g. a first derivative) and finally transformed
1830  back to the spatial domain with a sine transform (odd symmetric). For this to work
1831  properly the image must be shifted left or up by one pixel (depending on whether
1832  the x- or y-axis is odd symmetric) before the inverse transform can be applied.
1833  (see example below).
1834 
1835  The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
1836  where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
1837  whether the x- and y-axis is to be transformed using even or odd symmetry.
1838  The same functions can be used for both the forward and inverse transforms,
1839  only the normalization changes. For signal processing, the following
1840  normalization factors are most appropriate:
1841 
1842  \code
1843  forward inverse
1844  ------------------------------------------------------------
1845  X even, Y even 1.0 4.0 * (w-1) * (h-1)
1846  X even, Y odd -1.0 -4.0 * (w-1) * (h+1)
1847  X odd, Y even -1.0 -4.0 * (w+1) * (h-1)
1848  X odd, Y odd 1.0 4.0 * (w+1) * (h+1)
1849  \endcode
1850 
1851  where <tt>w</tt> and <tt>h</tt> denote the image width and height.
1852 
1853  <b> Declarations:</b>
1854 
1855  pass arguments explicitly:
1856  \code
1857  namespace vigra {
1858  template <class SrcTraverser, class SrcAccessor,
1859  class DestTraverser, class DestAccessor>
1860  void
1861  fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
1862  DestTraverser dul, DestAccessor dest, fftw_real norm);
1863 
1864  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
1865  }
1866  \endcode
1867 
1868 
1869  use argument objects in conjunction with \ref ArgumentObjectFactories :
1870  \code
1871  namespace vigra {
1872  template <class SrcTraverser, class SrcAccessor,
1873  class DestTraverser, class DestAccessor>
1874  void
1875  fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
1876  pair<DestTraverser, DestAccessor> dest, fftw_real norm);
1877 
1878  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
1879  }
1880  \endcode
1881 
1882  <b> Usage:</b>
1883 
1884  <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br>
1885  Namespace: vigra
1886 
1887  \code
1888  vigra::FImage spatial(width,height), fourier(width,height);
1889  ... // fill image with data
1890 
1891  // forward cosine transform == reflective boundary conditions
1892  fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
1893 
1894  // multiply with a first derivative of Gaussian in x-direction
1895  for(int y = 0; y < height; ++y)
1896  {
1897  for(int x = 1; x < width; ++x)
1898  {
1899  double dx = x * M_PI / (width - 1);
1900  double dy = y * M_PI / (height - 1);
1901  fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
1902  }
1903  fourier(width-1, y) = 0.0;
1904  }
1905 
1906  // inverse transform -- odd symmetry in x-direction, even in y,
1907  // due to symmetry of the filter
1908  fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
1909  (fftw_real)-4.0 * (width+1) * (height-1));
1910  \endcode
1911 */
1913 
1914 template <class SrcTraverser, class SrcAccessor,
1915  class DestTraverser, class DestAccessor>
1916 inline void
1917 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
1918  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
1919 {
1920  fourierTransformRealEE(src.first, src.second, src.third,
1921  dest.first, dest.second, norm);
1922 }
1923 
1924 template <class SrcTraverser, class SrcAccessor,
1925  class DestTraverser, class DestAccessor>
1926 inline void
1927 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
1928  DestTraverser dul, DestAccessor dest, fftw_real norm)
1929 {
1930  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
1931  norm, FFTW_REDFT00, FFTW_REDFT00);
1932 }
1933 
1934 template <class DestTraverser, class DestAccessor>
1935 inline void
1936 fourierTransformRealEE(
1940  DestTraverser dul, DestAccessor dest, fftw_real norm)
1941 {
1942  int w = slr.x - sul.x;
1943 
1944  // test for right memory layout (fftw expects a width*height fftw_real array)
1945  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
1946  fourierTransformRealImpl(sul, slr, dul, dest,
1947  norm, FFTW_REDFT00, FFTW_REDFT00);
1948  else
1949  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
1950  norm, FFTW_REDFT00, FFTW_REDFT00);
1951 }
1952 
1953 /********************************************************************/
1954 
1955 template <class SrcTraverser, class SrcAccessor,
1956  class DestTraverser, class DestAccessor>
1957 inline void
1958 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
1959  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
1960 {
1961  fourierTransformRealOE(src.first, src.second, src.third,
1962  dest.first, dest.second, norm);
1963 }
1964 
1965 template <class SrcTraverser, class SrcAccessor,
1966  class DestTraverser, class DestAccessor>
1967 inline void
1968 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
1969  DestTraverser dul, DestAccessor dest, fftw_real norm)
1970 {
1971  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
1972  norm, FFTW_RODFT00, FFTW_REDFT00);
1973 }
1974 
1975 template <class DestTraverser, class DestAccessor>
1976 inline void
1977 fourierTransformRealOE(
1981  DestTraverser dul, DestAccessor dest, fftw_real norm)
1982 {
1983  int w = slr.x - sul.x;
1984 
1985  // test for right memory layout (fftw expects a width*height fftw_real array)
1986  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
1987  fourierTransformRealImpl(sul, slr, dul, dest,
1988  norm, FFTW_RODFT00, FFTW_REDFT00);
1989  else
1990  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
1991  norm, FFTW_RODFT00, FFTW_REDFT00);
1992 }
1993 
1994 /********************************************************************/
1995 
1996 template <class SrcTraverser, class SrcAccessor,
1997  class DestTraverser, class DestAccessor>
1998 inline void
1999 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2000  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2001 {
2002  fourierTransformRealEO(src.first, src.second, src.third,
2003  dest.first, dest.second, norm);
2004 }
2005 
2006 template <class SrcTraverser, class SrcAccessor,
2007  class DestTraverser, class DestAccessor>
2008 inline void
2009 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2010  DestTraverser dul, DestAccessor dest, fftw_real norm)
2011 {
2012  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2013  norm, FFTW_REDFT00, FFTW_RODFT00);
2014 }
2015 
2016 template <class DestTraverser, class DestAccessor>
2017 inline void
2018 fourierTransformRealEO(
2022  DestTraverser dul, DestAccessor dest, fftw_real norm)
2023 {
2024  int w = slr.x - sul.x;
2025 
2026  // test for right memory layout (fftw expects a width*height fftw_real array)
2027  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2028  fourierTransformRealImpl(sul, slr, dul, dest,
2029  norm, FFTW_REDFT00, FFTW_RODFT00);
2030  else
2031  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2032  norm, FFTW_REDFT00, FFTW_RODFT00);
2033 }
2034 
2035 /********************************************************************/
2036 
2037 template <class SrcTraverser, class SrcAccessor,
2038  class DestTraverser, class DestAccessor>
2039 inline void
2040 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2041  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2042 {
2043  fourierTransformRealOO(src.first, src.second, src.third,
2044  dest.first, dest.second, norm);
2045 }
2046 
2047 template <class SrcTraverser, class SrcAccessor,
2048  class DestTraverser, class DestAccessor>
2049 inline void
2050 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2051  DestTraverser dul, DestAccessor dest, fftw_real norm)
2052 {
2053  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2054  norm, FFTW_RODFT00, FFTW_RODFT00);
2055 }
2056 
2057 template <class DestTraverser, class DestAccessor>
2058 inline void
2059 fourierTransformRealOO(
2063  DestTraverser dul, DestAccessor dest, fftw_real norm)
2064 {
2065  int w = slr.x - sul.x;
2066 
2067  // test for right memory layout (fftw expects a width*height fftw_real array)
2068  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2069  fourierTransformRealImpl(sul, slr, dul, dest,
2070  norm, FFTW_RODFT00, FFTW_RODFT00);
2071  else
2072  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2073  norm, FFTW_RODFT00, FFTW_RODFT00);
2074 }
2075 
2076 /*******************************************************************/
2077 
2078 template <class SrcTraverser, class SrcAccessor,
2079  class DestTraverser, class DestAccessor>
2080 void
2081 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2082  DestTraverser dul, DestAccessor dest,
2083  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2084 {
2085  FFTWRealImage workImage(slr - sul);
2086  copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2087  FFTWRealImage const & cworkImage = workImage;
2088  fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
2089  dul, dest, norm, kindx, kindy);
2090 }
2091 
2092 
2093 template <class DestTraverser, class DestAccessor>
2094 void
2095 fourierTransformRealImpl(
2098  DestTraverser dul, DestAccessor dest,
2099  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2100 {
2101  int w = slr.x - sul.x;
2102  int h = slr.y - sul.y;
2103  BasicImage<fftw_real> res(w, h);
2104 
2105  fftw_plan plan = fftw_plan_r2r_2d(h, w,
2106  (fftw_real *)&(*sul), (fftw_real *)res.begin(),
2107  kindy, kindx, FFTW_ESTIMATE);
2108  fftw_execute(plan);
2109  fftw_destroy_plan(plan);
2110 
2111  if(norm != 1.0)
2112  transformImage(srcImageRange(res), destIter(dul, dest),
2113  std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2114  else
2115  copyImage(srcImageRange(res), destIter(dul, dest));
2116 }
2117 
2118 
2119 //@}
2120 
2121 } // namespace vigra
2122 
2123 #endif // VIGRA_FFTW3_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.1 (Wed Mar 12 2014)