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

fftw3.hxx VIGRA

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 <complex>
42 #include "stdimage.hxx"
43 #include "copyimage.hxx"
44 #include "transformimage.hxx"
45 #include "combineimages.hxx"
46 #include "numerictraits.hxx"
47 #include "imagecontainer.hxx"
48 #include <fftw3.h>
49 
50 namespace vigra {
51 
52 typedef double fftw_real;
53 
54 template <class T>
55 struct FFTWReal;
56 
57 template <>
58 struct FFTWReal<fftw_complex>
59 {
60  typedef double type;
61 };
62 
63 template <>
64 struct FFTWReal<fftwf_complex>
65 {
66  typedef float type;
67 };
68 
69 template <>
70 struct FFTWReal<fftwl_complex>
71 {
72  typedef long double type;
73 };
74 
75 template <class T>
76 struct FFTWReal2Complex;
77 
78 template <>
79 struct FFTWReal2Complex<double>
80 {
81  typedef fftw_complex type;
82  typedef fftw_plan plan_type;
83 };
84 
85 template <>
86 struct FFTWReal2Complex<float>
87 {
88  typedef fftwf_complex type;
89  typedef fftwf_plan plan_type;
90 };
91 
92 template <>
93 struct FFTWReal2Complex<long double>
94 {
95  typedef fftwl_complex type;
96  typedef fftwl_plan plan_type;
97 };
98 
99 /********************************************************/
100 /* */
101 /* FFTWComplex */
102 /* */
103 /********************************************************/
104 
105 /** \brief Wrapper class for the FFTW complex types '<TT>fftw_complex</TT>'.
106 
107  This class encapsulates the low-level complex number types provided by the
108  <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a> library (i.e.
109  '<TT>fftw_complex</TT>', '<TT>fftwf_complex</TT>', '<TT>fftwl_complex</TT>').
110  In particular, it provides constructors, member functions and
111  \ref FFTWComplexOperators "arithmetic operators" that make FFTW complex numbers
112  compatible with <tt>std::complex</tt>. In addition, the class defines
113  transformations to polar coordinates and \ref FFTWComplexAccessors "accessors".
114 
115  FFTWComplex implements the concepts \ref AlgebraicField and
116  \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt>
117  and <tt>FFTWComplexImage</tt> are defined.
118 
119  <b>See also:</b>
120  <ul>
121  <li> \ref FFTWComplexTraits
122  <li> \ref FFTWComplexOperators
123  <li> \ref FFTWComplexAccessors
124  </ul>
125 
126  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
127  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
128  Namespace: vigra
129 */
130 template <class Real = double>
131 class FFTWComplex
132 {
133  public:
134  /** The wrapped complex type
135  */
136  typedef typename FFTWReal2Complex<Real>::type complex_type;
137 
138  /** The complex' component type, as defined in '<TT>fftw3.h</TT>'
139  */
140  typedef Real value_type;
141 
142  /** reference type (result of operator[])
143  */
144  typedef value_type & reference;
145 
146  /** const reference type (result of operator[] const)
147  */
148  typedef value_type const & const_reference;
149 
150  /** iterator type (result of begin() )
151  */
152  typedef value_type * iterator;
153 
154  /** const iterator type (result of begin() const)
155  */
156  typedef value_type const * const_iterator;
157 
158  /** The norm type (result of magnitude())
159  */
160  typedef value_type NormType;
161 
162  /** The squared norm type (result of squaredMagnitde())
163  */
164  typedef value_type SquaredNormType;
165 
166  /** Construct from real and imaginary part.
167  Default: 0.
168  */
169  FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0)
170  {
171  data_[0] = re;
172  data_[1] = im;
173  }
174 
175  /** Copy constructor.
176  */
177  FFTWComplex(FFTWComplex const & o)
178  {
179  data_[0] = o.data_[0];
180  data_[1] = o.data_[1];
181  }
182 
183  /** Copy constructor.
184  */
185  template <class U>
186  FFTWComplex(FFTWComplex<U> const & o)
187  {
188  data_[0] = (Real)o.real();
189  data_[1] = (Real)o.imag();
190  }
191 
192  /** Construct from plain <TT>fftw_complex</TT>.
193  */
194  FFTWComplex(fftw_complex const & o)
195  {
196  data_[0] = (Real)o[0];
197  data_[1] = (Real)o[1];
198  }
199 
200  /** Construct from plain <TT>fftwf_complex</TT>.
201  */
202  FFTWComplex(fftwf_complex const & o)
203  {
204  data_[0] = (Real)o[0];
205  data_[1] = (Real)o[1];
206  }
207 
208  /** Construct from plain <TT>fftwl_complex</TT>.
209  */
210  FFTWComplex(fftwl_complex const & o)
211  {
212  data_[0] = (Real)o[0];
213  data_[1] = (Real)o[1];
214  }
215 
216  /** Construct from std::complex.
217  */
218  template <class T>
219  FFTWComplex(std::complex<T> const & o)
220  {
221  data_[0] = (Real)o.real();
222  data_[1] = (Real)o.imag();
223  }
224 
225  /** Construct from TinyVector.
226  */
227  template <class T>
228  FFTWComplex(TinyVector<T, 2> const & o)
229  {
230  data_[0] = (Real)o[0];
231  data_[1] = (Real)o[1];
232  }
233 
234  /** Assignment.
235  */
236  FFTWComplex& operator=(FFTWComplex const & o)
237  {
238  data_[0] = o.data_[0];
239  data_[1] = o.data_[1];
240  return *this;
241  }
242 
243  /** Assignment.
244  */
245  template <class U>
246  FFTWComplex& operator=(FFTWComplex<U> const & o)
247  {
248  data_[0] = (Real)o.real();
249  data_[1] = (Real)o.imag();
250  return *this;
251  }
252 
253  /** Assignment.
254  */
255  FFTWComplex& operator=(fftw_complex const & o)
256  {
257  data_[0] = (Real)o[0];
258  data_[1] = (Real)o[1];
259  return *this;
260  }
261 
262  /** Assignment.
263  */
264  FFTWComplex& operator=(fftwf_complex const & o)
265  {
266  data_[0] = (Real)o[0];
267  data_[1] = (Real)o[1];
268  return *this;
269  }
270 
271  /** Assignment.
272  */
273  FFTWComplex& operator=(fftwl_complex const & o)
274  {
275  data_[0] = (Real)o[0];
276  data_[1] = (Real)o[1];
277  return *this;
278  }
279 
280  /** Assignment.
281  */
282  FFTWComplex& operator=(double o)
283  {
284  data_[0] = (Real)o;
285  data_[1] = 0.0;
286  return *this;
287  }
288 
289  /** Assignment.
290  */
291  FFTWComplex& operator=(float o)
292  {
293  data_[0] = (Real)o;
294  data_[1] = 0.0;
295  return *this;
296  }
297 
298  /** Assignment.
299  */
300  FFTWComplex& operator=(long double o)
301  {
302  data_[0] = (Real)o;
303  data_[1] = 0.0;
304  return *this;
305  }
306 
307  /** Assignment.
308  */
309  template <class T>
310  FFTWComplex& operator=(TinyVector<T, 2> const & o)
311  {
312  data_[0] = (Real)o[0];
313  data_[1] = (Real)o[1];
314  return *this;
315  }
316 
317  /** Assignment.
318  */
319  template <class T>
320  FFTWComplex& operator=(std::complex<T> const & o)
321  {
322  data_[0] = (Real)o.real();
323  data_[1] = (Real)o.imag();
324  return *this;
325  }
326 
327  reference re()
328  { return data_[0]; }
329 
330  const_reference re() const
331  { return data_[0]; }
332 
333  reference real()
334  { return data_[0]; }
335 
336  const_reference real() const
337  { return data_[0]; }
338 
339  reference im()
340  { return data_[1]; }
341 
342  const_reference im() const
343  { return data_[1]; }
344 
345  reference imag()
346  { return data_[1]; }
347 
348  const_reference imag() const
349  { return data_[1]; }
350 
351  /** Unary negation.
352  */
353  FFTWComplex operator-() const
354  { return FFTWComplex(-data_[0], -data_[1]); }
355 
356  /** Squared magnitude x*conj(x)
357  */
358  SquaredNormType squaredMagnitude() const
359  { return data_[0]*data_[0]+data_[1]*data_[1]; }
360 
361  /** Magnitude (length of radius vector).
362  */
363  NormType magnitude() const
364  { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
365 
366  /** Phase angle.
367  */
368  value_type phase() const
369  { return VIGRA_CSTD::atan2(data_[1], data_[0]); }
370 
371  /** Access components as if number were a vector.
372  */
373  reference operator[](int i)
374  { return data_[i]; }
375 
376  /** Read components as if number were a vector.
377  */
378  const_reference operator[](int i) const
379  { return data_[i]; }
380 
381  /** Length of complex number (always 2).
382  */
383  int size() const
384  { return 2; }
385 
386  iterator begin()
387  { return data_; }
388 
389  iterator end()
390  { return data_ + 2; }
391 
392  const_iterator begin() const
393  { return data_; }
394 
395  const_iterator end() const
396  { return data_ + 2; }
397 
398  private:
399  complex_type data_;
400 };
401 
402 /********************************************************/
403 /* */
404 /* FFTWComplexTraits */
405 /* */
406 /********************************************************/
407 
408 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex
409 
410  The numeric and promote traits for fftw_complex and FFTWComplex follow
411  the general specifications for \ref NumericPromotionTraits and
412  \ref AlgebraicField. They are explicitly specialized for the types
413  involved:
414 
415  \code
416 
417  template<>
418  struct NumericTraits<fftw_complex>
419  {
420  typedef fftw_complex Promote;
421  typedef fftw_complex RealPromote;
422  typedef fftw_complex ComplexPromote;
423  typedef fftw_real ValueType;
424 
425  typedef VigraFalseType isIntegral;
426  typedef VigraFalseType isScalar;
427  typedef VigraFalseType isOrdered;
428  typedef VigraTrueType isComplex;
429 
430  // etc.
431  };
432 
433  template<class Real>
434  struct NumericTraits<FFTWComplex<Real> >
435  {
436  typedef FFTWComplex<Real> Promote;
437  typedef FFTWComplex<Real> RealPromote;
438  typedef FFTWComplex<Real> ComplexPromote;
439  typedef Real ValueType;
440 
441  typedef VigraFalseType isIntegral;
442  typedef VigraFalseType isScalar;
443  typedef VigraFalseType isOrdered;
444  typedef VigraTrueType isComplex;
445 
446  // etc.
447  };
448 
449  template<>
450  struct NormTraits<fftw_complex>
451  {
452  typedef fftw_complex Type;
453  typedef fftw_real SquaredNormType;
454  typedef fftw_real NormType;
455  };
456 
457  template<class Real>
458  struct NormTraits<FFTWComplex>
459  {
460  typedef FFTWComplex Type;
461  typedef fftw_real SquaredNormType;
462  typedef fftw_real NormType;
463  };
464 
465  template <>
466  struct PromoteTraits<fftw_complex, fftw_complex>
467  {
468  typedef fftw_complex Promote;
469  };
470 
471  template <>
472  struct PromoteTraits<fftw_complex, double>
473  {
474  typedef fftw_complex Promote;
475  };
476 
477  template <>
478  struct PromoteTraits<double, fftw_complex>
479  {
480  typedef fftw_complex Promote;
481  };
482 
483  template <class Real>
484  struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
485  {
486  typedef FFTWComplex<Real> Promote;
487  };
488 
489  template <class Real>
490  struct PromoteTraits<FFTWComplex<Real>, double>
491  {
492  typedef FFTWComplex<Real> Promote;
493  };
494 
495  template <class Real>
496  struct PromoteTraits<double, FFTWComplex<Real> >
497  {
498  typedef FFTWComplex<Real> Promote;
499  };
500  \endcode
501 
502  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
503  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
504  Namespace: vigra
505 
506 */
507 template<>
508 struct NumericTraits<fftw_complex>
509 {
510  typedef fftw_complex Type;
511  typedef fftw_complex Promote;
512  typedef fftw_complex RealPromote;
513  typedef fftw_complex ComplexPromote;
514  typedef fftw_real ValueType;
515 
516  typedef VigraFalseType isIntegral;
517  typedef VigraFalseType isScalar;
518  typedef NumericTraits<fftw_real>::isSigned isSigned;
519  typedef VigraFalseType isOrdered;
520  typedef VigraTrueType isComplex;
521 
522  static FFTWComplex<> zero() { return FFTWComplex<>(0.0, 0.0); }
523  static FFTWComplex<> one() { return FFTWComplex<>(1.0, 0.0); }
524  static FFTWComplex<> nonZero() { return one(); }
525 
526  static const Promote & toPromote(const Type & v) { return v; }
527  static const RealPromote & toRealPromote(const Type & v) { return v; }
528  static const Type & fromPromote(const Promote & v) { return v; }
529  static const Type & fromRealPromote(const RealPromote & v) { return v; }
530 };
531 
532 template<class Real>
533 struct NumericTraits<FFTWComplex<Real> >
534 {
535  typedef FFTWComplex<Real> Type;
536  typedef FFTWComplex<Real> Promote;
537  typedef FFTWComplex<Real> RealPromote;
538  typedef FFTWComplex<Real> ComplexPromote;
539  typedef typename Type::value_type ValueType;
540 
541  typedef VigraFalseType isIntegral;
542  typedef VigraFalseType isScalar;
543  typedef typename NumericTraits<ValueType>::isSigned isSigned;
544  typedef VigraFalseType isOrdered;
545  typedef VigraTrueType isComplex;
546 
547  static Type zero() { return Type(0.0, 0.0); }
548  static Type one() { return Type(1.0, 0.0); }
549  static Type nonZero() { return one(); }
550 
551  static const Promote & toPromote(const Type & v) { return v; }
552  static const RealPromote & toRealPromote(const Type & v) { return v; }
553  static const Type & fromPromote(const Promote & v) { return v; }
554  static const Type & fromRealPromote(const RealPromote & v) { return v; }
555 };
556 
557 template<>
558 struct NormTraits<fftw_complex>
559 {
560  typedef fftw_complex Type;
561  typedef fftw_real SquaredNormType;
562  typedef fftw_real NormType;
563 };
564 
565 template<class Real>
566 struct NormTraits<FFTWComplex<Real> >
567 {
568  typedef FFTWComplex<Real> Type;
569  typedef typename Type::SquaredNormType SquaredNormType;
570  typedef typename Type::NormType NormType;
571 };
572 
573 template <>
574 struct PromoteTraits<fftw_complex, fftw_complex>
575 {
576  typedef fftw_complex Promote;
577 };
578 
579 template <>
580 struct PromoteTraits<fftw_complex, double>
581 {
582  typedef fftw_complex Promote;
583 };
584 
585 template <>
586 struct PromoteTraits<double, fftw_complex>
587 {
588  typedef fftw_complex Promote;
589 };
590 
591 template <class Real>
592 struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
593 {
594  typedef FFTWComplex<Real> Promote;
595 };
596 
597 template <class Real>
598 struct PromoteTraits<FFTWComplex<Real>, double>
599 {
600  typedef FFTWComplex<Real> Promote;
601 };
602 
603 template <class Real>
604 struct PromoteTraits<double, FFTWComplex<Real> >
605 {
606  typedef FFTWComplex<Real> Promote;
607 };
608 
609 template<class T>
610 struct CanSkipInitialization<std::complex<T> >
611 {
612  typedef typename CanSkipInitialization<T>::type type;
613  static const bool value = type::asBool;
614 };
615 
616 template<class Real>
617 struct CanSkipInitialization<FFTWComplex<Real> >
618 {
619  typedef typename CanSkipInitialization<Real>::type type;
620  static const bool value = type::asBool;
621 };
622 
623 namespace multi_math {
624 
625 template <class ARG>
626 struct MultiMathOperand;
627 
628 template <class Real>
629 struct MultiMathOperand<FFTWComplex<Real> >
630 {
631  typedef MultiMathOperand<FFTWComplex<Real> > AllowOverload;
632  typedef FFTWComplex<Real> result_type;
633 
634  static const int ndim = 0;
635 
636  MultiMathOperand(FFTWComplex<Real> const & v)
637  : v_(v)
638  {}
639 
640  template <class SHAPE>
641  bool checkShape(SHAPE const &) const
642  {
643  return true;
644  }
645 
646  template <class SHAPE>
647  FFTWComplex<Real> const & operator[](SHAPE const &) const
648  {
649  return v_;
650  }
651 
652  void inc(unsigned int /*LEVEL*/) const
653  {}
654 
655  void reset(unsigned int /*LEVEL*/) const
656  {}
657 
658  FFTWComplex<Real> const & operator*() const
659  {
660  return v_;
661  }
662 
663  FFTWComplex<Real> v_;
664 };
665 
666 } // namespace multi_math
667 
668 template<class Ty>
669 class FFTWAllocator
670 {
671  public:
672  typedef std::size_t size_type;
673  typedef std::ptrdiff_t difference_type;
674  typedef Ty *pointer;
675  typedef const Ty *const_pointer;
676  typedef Ty& reference;
677  typedef const Ty& const_reference;
678  typedef Ty value_type;
679 
680  pointer address(reference val) const
681  { return &val; }
682 
683  const_pointer address(const_reference val) const
684  { return &val; }
685 
686  template<class Other>
687  struct rebind
688  {
689  typedef FFTWAllocator<Other> other;
690  };
691 
692  FFTWAllocator() throw()
693  {}
694 
695  template<class Other>
696  FFTWAllocator(const FFTWAllocator<Other>& /*right*/) throw()
697  {}
698 
699  template<class Other>
700  FFTWAllocator& operator=(const FFTWAllocator<Other>& /*right*/)
701  {
702  return *this;
703  }
704 
705  pointer allocate(size_type count, void * = 0)
706  {
707  return (pointer)fftw_malloc(count * sizeof(Ty));
708  }
709 
710  void deallocate(pointer ptr, size_type /*count*/)
711  {
712  fftw_free(ptr);
713  }
714 
715  void construct(pointer ptr, const Ty& val)
716  {
717  new(ptr) Ty(val);
718 
719  }
720 
721  void destroy(pointer ptr)
722  {
723  ptr->~Ty();
724  }
725 
726  size_type max_size() const throw()
727  {
728  return NumericTraits<std::ptrdiff_t>::max() / sizeof(Ty);
729  }
730 };
731 
732 } // namespace vigra
733 
734 namespace std {
735 
736 template<class Real>
737 class allocator<vigra::FFTWComplex<Real> >
738 {
739  public:
740  typedef vigra::FFTWComplex<Real> value_type;
741  typedef size_t size_type;
742  typedef std::ptrdiff_t difference_type;
743  typedef value_type *pointer;
744  typedef const value_type *const_pointer;
745  typedef value_type& reference;
746  typedef const value_type& const_reference;
747 
748  pointer address(reference val) const
749  { return &val; }
750 
751  const_pointer address(const_reference val) const
752  { return &val; }
753 
754  template<class Other>
755  struct rebind
756  {
757  typedef allocator<Other> other;
758  };
759 
760  allocator() throw()
761  {}
762 
763  template<class Other>
764  allocator(const allocator<Other>& /*right*/) throw()
765  {}
766 
767  template<class Other>
768  allocator& operator=(const allocator<Other>& /*right*/)
769  {
770  return *this;
771  }
772 
773  pointer allocate(size_type count, void * = 0)
774  {
775  return (pointer)fftw_malloc(count * sizeof(value_type));
776  }
777 
778  void deallocate(pointer ptr, size_type /*count*/)
779  {
780  fftw_free(ptr);
781  }
782 
783  void construct(pointer ptr, const value_type& val)
784  {
785  new(ptr) value_type(val);
786 
787  }
788 
789  void destroy(pointer ptr)
790  {
791  ptr->~value_type();
792  }
793 
794  size_type max_size() const throw()
795  {
796  return vigra::NumericTraits<std::ptrdiff_t>::max() / sizeof(value_type);
797  }
798 };
799 
800 } // namespace std
801 
802 namespace vigra {
803 
804 /********************************************************/
805 /* */
806 /* FFTWComplex Operations */
807 /* */
808 /********************************************************/
809 
810 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex
811 
812  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
813  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
814 
815  These functions fulfill the requirements of an Algebraic Field.
816  Return types are determined according to \ref FFTWComplexTraits.
817 
818  Namespace: vigra
819  <p>
820 
821  */
822 //@{
823  /// equal
824 template <class R>
825 inline bool operator ==(FFTWComplex<R> const &a, const FFTWComplex<R> &b) {
826  return a.re() == b.re() && a.im() == b.im();
827 }
828 
829 template <class R>
830 inline bool operator ==(FFTWComplex<R> const &a, double b) {
831  return a.re() == b && a.im() == 0.0;
832 }
833 
834 template <class R>
835 inline bool operator ==(double a, const FFTWComplex<R> &b) {
836  return a == b.re() && 0.0 == b.im();
837 }
838 
839  /// not equal
840 template <class R>
841 inline bool operator !=(FFTWComplex<R> const &a, const FFTWComplex<R> &b) {
842  return a.re() != b.re() || a.im() != b.im();
843 }
844 
845  /// not equal
846 template <class R>
847 inline bool operator !=(FFTWComplex<R> const &a, double b) {
848  return a.re() != b || a.im() != 0.0;
849 }
850 
851  /// not equal
852 template <class R>
853 inline bool operator !=(double a, const FFTWComplex<R> &b) {
854  return a != b.re() || 0.0 != b.im();
855 }
856 
857  /// add-assignment
858 template <class R>
859 inline FFTWComplex<R> & operator +=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
860  a.re() += b.re();
861  a.im() += b.im();
862  return a;
863 }
864 
865  /// subtract-assignment
866 template <class R>
867 inline FFTWComplex<R> & operator -=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
868  a.re() -= b.re();
869  a.im() -= b.im();
870  return a;
871 }
872 
873  /// multiply-assignment
874 template <class R>
875 inline FFTWComplex<R> & operator *=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
876  typename FFTWComplex<R>::value_type t = a.re()*b.re()-a.im()*b.im();
877  a.im() = a.re()*b.im()+a.im()*b.re();
878  a.re() = t;
879  return a;
880 }
881 
882  /// divide-assignment
883 template <class R>
884 inline FFTWComplex<R> & operator /=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
885  typename FFTWComplex<R>::value_type sm = b.squaredMagnitude();
886  typename FFTWComplex<R>::value_type t = (a.re()*b.re()+a.im()*b.im())/sm;
887  a.im() = (b.re()*a.im()-a.re()*b.im())/sm;
888  a.re() = t;
889  return a;
890 }
891 
892  /// add-assignment with scalar double
893 template <class R>
894 inline FFTWComplex<R> & operator +=(FFTWComplex<R> &a, double b) {
895  a.re() += (R)b;
896  return a;
897 }
898 
899  /// subtract-assignment with scalar double
900 template <class R>
901 inline FFTWComplex<R> & operator -=(FFTWComplex<R> &a, double b) {
902  a.re() -= (R)b;
903  return a;
904 }
905 
906  /// multiply-assignment with scalar double
907 template <class R>
908 inline FFTWComplex<R> & operator *=(FFTWComplex<R> &a, double b) {
909  a.re() *= (R)b;
910  a.im() *= (R)b;
911  return a;
912 }
913 
914  /// divide-assignment with scalar double
915 template <class R>
916 inline FFTWComplex<R> & operator /=(FFTWComplex<R> &a, double b) {
917  a.re() /= (R)b;
918  a.im() /= (R)b;
919  return a;
920 }
921 
922  /// addition
923 template <class R>
924 inline FFTWComplex<R> operator +(FFTWComplex<R> a, const FFTWComplex<R> &b) {
925  a += b;
926  return a;
927 }
928 
929  /// right addition with scalar double
930 template <class R>
931 inline FFTWComplex<R> operator +(FFTWComplex<R> a, double b) {
932  a += b;
933  return a;
934 }
935 
936  /// left addition with scalar double
937 template <class R>
938 inline FFTWComplex<R> operator +(double a, FFTWComplex<R> b) {
939  b += a;
940  return b;
941 }
942 
943  /// subtraction
944 template <class R>
945 inline FFTWComplex<R> operator -(FFTWComplex<R> a, const FFTWComplex<R> &b) {
946  a -= b;
947  return a;
948 }
949 
950  /// right subtraction with scalar double
951 template <class R>
952 inline FFTWComplex<R> operator -(FFTWComplex<R> a, double b) {
953  a -= b;
954  return a;
955 }
956 
957  /// left subtraction with scalar double
958 template <class R>
959 inline FFTWComplex<R> operator -(double a, FFTWComplex<R> const & b) {
960  return (-b) += a;
961 }
962 
963  /// multiplication
964 template <class R>
965 inline FFTWComplex<R> operator *(FFTWComplex<R> a, const FFTWComplex<R> &b) {
966  a *= b;
967  return a;
968 }
969 
970  /// right multiplication with scalar double
971 template <class R>
972 inline FFTWComplex<R> operator *(FFTWComplex<R> a, double b) {
973  a *= b;
974  return a;
975 }
976 
977  /// left multiplication with scalar double
978 template <class R>
979 inline FFTWComplex<R> operator *(double a, FFTWComplex<R> b) {
980  b *= a;
981  return b;
982 }
983 
984  /// division
985 template <class R>
986 inline FFTWComplex<R> operator /(FFTWComplex<R> a, const FFTWComplex<R> &b) {
987  a /= b;
988  return a;
989 }
990 
991  /// right division with scalar double
992 template <class R>
993 inline FFTWComplex<R> operator /(FFTWComplex<R> a, double b) {
994  a /= b;
995  return a;
996 }
997 
998 using VIGRA_CSTD::abs;
999 
1000  /// absolute value (= magnitude)
1001 template <class R>
1002 inline typename FFTWComplex<R>::NormType abs(const FFTWComplex<R> &a)
1003 {
1004  return a.magnitude();
1005 }
1006 
1007  /// phase
1008 template <class R>
1009 inline R arg(const FFTWComplex<R> &a)
1010 {
1011  return a.phase();
1012 }
1013 
1014  /// real part
1015 template <class R>
1016 inline R real(const FFTWComplex<R> &a)
1017 {
1018  return a.real();
1019 }
1020 
1021  /// imaginary part
1022 template <class R>
1023 inline R imag(const FFTWComplex<R> &a)
1024 {
1025  return a.imag();
1026 }
1027 
1028  /// complex conjugate
1029 template <class R>
1030 inline FFTWComplex<R> conj(const FFTWComplex<R> &a)
1031 {
1032  return FFTWComplex<R>(a.re(), -a.im());
1033 }
1034 
1035  /// norm (= magnitude)
1036 template <class R>
1037 inline typename FFTWComplex<R>::NormType norm(const FFTWComplex<R> &a)
1038 {
1039  return a.magnitude();
1040 }
1041 
1042  /// squared norm (= squared magnitude)
1043 template <class R>
1044 inline typename FFTWComplex<R>::SquaredNormType squaredNorm(const FFTWComplex<R> &a)
1045 {
1046  return a.squaredMagnitude();
1047 }
1048 
1049 #define VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(fct) \
1050 template <class R> \
1051 inline FFTWComplex<R> fct(const FFTWComplex<R> &a) \
1052 { \
1053  return std::fct(reinterpret_cast<std::complex<R> const &>(a)); \
1054 }
1055 
1056 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cos)
1057 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cosh)
1058 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(exp)
1059 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(log)
1060 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(log10)
1061 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sin)
1062 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sinh)
1063 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sqrt)
1064 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tan)
1065 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tanh)
1066 
1067 #undef VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION
1068 
1069 template <class R>
1070 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, int e)
1071 {
1072  return std::pow(reinterpret_cast<std::complex<R> const &>(a), e);
1073 }
1074 
1075 template <class R>
1076 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, R const & e)
1077 {
1078  return std::pow(reinterpret_cast<std::complex<R> const &>(a), e);
1079 }
1080 
1081 template <class R>
1082 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, const FFTWComplex<R> & e)
1083 {
1084  return std::pow(reinterpret_cast<std::complex<R> const &>(a),
1085  reinterpret_cast<std::complex<R> const &>(e));
1086 }
1087 
1088 template <class R>
1089 inline FFTWComplex<R> pow(R const & a, const FFTWComplex<R> &e)
1090 {
1091  return std::pow(a, reinterpret_cast<std::complex<R> const &>(e));
1092 }
1093 
1094 //@}
1095 
1096 } // namespace vigra
1097 
1098 namespace std {
1099 
1100 template <class Real>
1101 ostream & operator<<(ostream & s, vigra::FFTWComplex<Real> const & v)
1102 {
1103  s << std::complex<Real>(v.re(), v.im());
1104  return s;
1105 }
1106 
1107 } // namespace std
1108 
1109 namespace vigra {
1110 
1111 /** \addtogroup StandardImageTypes
1112 */
1113 //@{
1114 
1115 /********************************************************/
1116 /* */
1117 /* FFTWRealImage */
1118 /* */
1119 /********************************************************/
1120 
1121  /** Float (<tt>fftw_real</tt>) image.
1122 
1123  The type <tt>fftw_real</tt> is defined as <tt>double</tt> (in FFTW 2 it used to be
1124  either <tt>float</tt> or <tt>double</tt>, as specified during compilation of FFTW).
1125  FFTWRealImage uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
1126  their const counterparts to access the data.
1127 
1128  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1129  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1130  Namespace: vigra
1131  */
1132 typedef BasicImage<fftw_real> FFTWRealImage;
1133 
1134 /********************************************************/
1135 /* */
1136 /* FFTWComplexImage */
1137 /* */
1138 /********************************************************/
1139 
1140 template<class R>
1141 struct IteratorTraits<
1142  BasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> >
1143 {
1144  typedef FFTWComplex<R> Type;
1145  typedef BasicImageIterator<Type, Type **> Iterator;
1146  typedef Iterator iterator;
1147  typedef BasicImageIterator<Type, Type **> mutable_iterator;
1148  typedef ConstBasicImageIterator<Type, Type **> const_iterator;
1149  typedef typename iterator::iterator_category iterator_category;
1150  typedef typename iterator::value_type value_type;
1151  typedef typename iterator::reference reference;
1152  typedef typename iterator::index_reference index_reference;
1153  typedef typename iterator::pointer pointer;
1154  typedef typename iterator::difference_type difference_type;
1155  typedef typename iterator::row_iterator row_iterator;
1156  typedef typename iterator::column_iterator column_iterator;
1157  typedef VectorAccessor<Type> default_accessor;
1158  typedef VectorAccessor<Type> DefaultAccessor;
1159  typedef VigraTrueType hasConstantStrides;
1160 };
1161 
1162 template<class R>
1163 struct IteratorTraits<
1164  ConstBasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> >
1165 {
1166  typedef FFTWComplex<R> Type;
1167  typedef ConstBasicImageIterator<Type, Type **> Iterator;
1168  typedef Iterator iterator;
1169  typedef BasicImageIterator<Type, Type **> mutable_iterator;
1170  typedef ConstBasicImageIterator<Type, Type **> const_iterator;
1171  typedef typename iterator::iterator_category iterator_category;
1172  typedef typename iterator::value_type value_type;
1173  typedef typename iterator::reference reference;
1174  typedef typename iterator::index_reference index_reference;
1175  typedef typename iterator::pointer pointer;
1176  typedef typename iterator::difference_type difference_type;
1177  typedef typename iterator::row_iterator row_iterator;
1178  typedef typename iterator::column_iterator column_iterator;
1179  typedef VectorAccessor<Type> default_accessor;
1180  typedef VectorAccessor<Type> DefaultAccessor;
1181  typedef VigraTrueType hasConstantStrides;
1182 };
1183 
1184  /** Complex (FFTWComplex) image.
1185  It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
1186  their const counterparts to access the data.
1187 
1188  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1189  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1190  Namespace: vigra
1191  */
1192 typedef BasicImage<FFTWComplex<> > FFTWComplexImage;
1193 
1194 //@}
1195 
1196 /********************************************************/
1197 /* */
1198 /* FFTWComplex-Accessors */
1199 /* */
1200 /********************************************************/
1201 
1202 /** \addtogroup DataAccessors
1203 */
1204 //@{
1205 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex
1206 
1207  Encapsulate access to pixels of type FFTWComplex
1208 */
1209 //@{
1210  /** Encapsulate access to the the real part of a complex number.
1211 
1212  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1213  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1214  Namespace: vigra
1215  */
1216 template <class Real = double>
1217 class FFTWRealAccessor
1218 {
1219  public:
1220 
1221  /// The accessor's value type.
1222  typedef Real value_type;
1223 
1224  /// Read real part at iterator position.
1225  template <class ITERATOR>
1226  value_type operator()(ITERATOR const & i) const {
1227  return (*i).re();
1228  }
1229 
1230  /// Read real part at offset from iterator position.
1231  template <class ITERATOR, class DIFFERENCE>
1232  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1233  return i[d].re();
1234  }
1235 
1236  /// Write real part at iterator position from a scalar.
1237  template <class ITERATOR>
1238  void set(value_type const & v, ITERATOR const & i) const {
1239  (*i).re()= v;
1240  }
1241 
1242  /// Write real part at offset from iterator position from a scalar.
1243  template <class ITERATOR, class DIFFERENCE>
1244  void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
1245  i[d].re()= v;
1246  }
1247 
1248  /// Write real part at iterator position into a scalar.
1249  template <class R, class ITERATOR>
1250  void set(FFTWComplex<R> const & v, ITERATOR const & i) const {
1251  *i = v.re();
1252  }
1253 
1254  /// Write real part at offset from iterator position into a scalar.
1255  template <class R, class ITERATOR, class DIFFERENCE>
1256  void set(FFTWComplex<R> const & v, ITERATOR const & i, DIFFERENCE d) const {
1257  i[d] = v.re();
1258  }
1259 };
1260 
1261  /** Encapsulate access to the the imaginary part of a complex number.
1262 
1263  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1264  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1265  Namespace: vigra
1266  */
1267 template <class Real = double>
1268 class FFTWImaginaryAccessor
1269 {
1270  public:
1271  /// The accessor's value type.
1272  typedef Real value_type;
1273 
1274  /// Read imaginary part at iterator position.
1275  template <class ITERATOR>
1276  value_type operator()(ITERATOR const & i) const {
1277  return (*i).im();
1278  }
1279 
1280  /// Read imaginary part at offset from iterator position.
1281  template <class ITERATOR, class DIFFERENCE>
1282  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1283  return i[d].im();
1284  }
1285 
1286  /// Write imaginary part at iterator position from a scalar.
1287  template <class ITERATOR>
1288  void set(value_type const & v, ITERATOR const & i) const {
1289  (*i).im()= v;
1290  }
1291 
1292  /// Write imaginary part at offset from iterator position from a scalar.
1293  template <class ITERATOR, class DIFFERENCE>
1294  void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
1295  i[d].im()= v;
1296  }
1297 
1298  /// Write imaginary part at iterator position into a scalar.
1299  template <class R, class ITERATOR>
1300  void set(FFTWComplex<R> const & v, ITERATOR const & i) const {
1301  *i = v.im();
1302  }
1303 
1304  /// Write imaginary part at offset from iterator position into a scalar.
1305  template <class R, class ITERATOR, class DIFFERENCE>
1306  void set(FFTWComplex<R> const & v, ITERATOR const & i, DIFFERENCE d) const {
1307  i[d] = v.im();
1308  }
1309 };
1310 
1311  /** Write a real number into a complex one. The imaginary part is set
1312  to 0.
1313 
1314  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1315  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1316  Namespace: vigra
1317  */
1318 template <class Real = double>
1319 class FFTWWriteRealAccessor
1320 : public FFTWRealAccessor<Real>
1321 {
1322  public:
1323  /// The accessor's value type.
1324  typedef Real value_type;
1325 
1326  /** Write real number at iterator position. Set imaginary part
1327  to 0.
1328  */
1329  template <class ITERATOR>
1330  void set(value_type const & v, ITERATOR const & i) const {
1331  (*i).re()= v;
1332  (*i).im()= 0;
1333  }
1334 
1335  /** Write real number at offset from iterator position. Set imaginary part
1336  to 0.
1337  */
1338  template <class ITERATOR, class DIFFERENCE>
1339  void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
1340  i[d].re()= v;
1341  i[d].im()= 0;
1342  }
1343 };
1344 
1345  /** Calculate squared magnitude of complex number on the fly.
1346 
1347  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1348  Namespace: vigra
1349  */
1350 template <class Real = double>
1351 class FFTWSquaredMagnitudeAccessor
1352 {
1353  public:
1354  /// The accessor's value type.
1355  typedef Real value_type;
1356 
1357  /// Read squared magnitude at iterator position.
1358  template <class ITERATOR>
1359  value_type operator()(ITERATOR const & i) const {
1360  return (*i).squaredMagnitude();
1361  }
1362 
1363  /// Read squared magnitude at offset from iterator position.
1364  template <class ITERATOR, class DIFFERENCE>
1365  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1366  return (i[d]).squaredMagnitude();
1367  }
1368 };
1369 
1370  /** Calculate magnitude of complex number on the fly.
1371 
1372  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1373  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1374  Namespace: vigra
1375  */
1376 template <class Real = double>
1377 class FFTWMagnitudeAccessor
1378 {
1379  public:
1380  /// The accessor's value type.
1381  typedef Real value_type;
1382 
1383  /// Read magnitude at iterator position.
1384  template <class ITERATOR>
1385  value_type operator()(ITERATOR const & i) const {
1386  return (*i).magnitude();
1387  }
1388 
1389  /// Read magnitude at offset from iterator position.
1390  template <class ITERATOR, class DIFFERENCE>
1391  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1392  return (i[d]).magnitude();
1393  }
1394 };
1395 
1396  /** Calculate natural logarithm of magnitude of complex number on the fly.
1397 
1398  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1399  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1400  Namespace: vigra
1401  */
1402 template <class Real = double>
1403 class FFTWLogMagnitudeAccessor
1404 {
1405  public:
1406  /// The accessor's value type.
1407  typedef Real value_type;
1408 
1409  /// Read natural log of magnitude at iterator position.
1410  template <class ITERATOR>
1411  value_type operator()(ITERATOR const & i) const {
1412  return std::log((*i).magnitude() + 1);
1413  }
1414 
1415  /// Read natural log of magnitude at offset from iterator position.
1416  template <class ITERATOR, class DIFFERENCE>
1417  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1418  return std::log((i[d]).magnitude() + 1);
1419  }
1420 };
1421 
1422  /** Calculate phase of complex number on the fly.
1423 
1424  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
1425  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
1426  Namespace: vigra
1427  */
1428 template <class Real = double>
1429 class FFTWPhaseAccessor
1430 {
1431  public:
1432  /// The accessor's value type.
1433  typedef Real value_type;
1434 
1435  /// Read phase at iterator position.
1436  template <class ITERATOR>
1437  value_type operator()(ITERATOR const & i) const {
1438  return (*i).phase();
1439  }
1440 
1441  /// Read phase at offset from iterator position.
1442  template <class ITERATOR, class DIFFERENCE>
1443  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1444  return (i[d]).phase();
1445  }
1446 };
1447 
1448 //@}
1449 //@}
1450 
1451 /********************************************************/
1452 /* */
1453 /* Fourier Transform */
1454 /* */
1455 /********************************************************/
1456 
1457 /* \addtogroup FourierTransform Fast Fourier Transform
1458 
1459  This documentation describes the VIGRA interface to FFTW version 3. The interface
1460  to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated.
1461 
1462  VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
1463  Transform</a> package to perform Fourier transformations. VIGRA
1464  provides a wrapper for FFTW's complex number type (FFTWComplex),
1465  but FFTW's functions are used verbatim. If the image is stored as
1466  a FFTWComplexImage, the simplest call to an FFT function is like this:
1467 
1468  \code
1469  vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
1470  ... // fill image with data
1471 
1472  // create a plan with estimated performance optimization
1473  fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
1474  (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
1475  FFTW_FORWARD, FFTW_ESTIMATE );
1476  // calculate FFT (this can be repeated as often as needed,
1477  // with fresh data written into the source array)
1478  fftw_execute(forwardPlan);
1479 
1480  // release the plan memory
1481  fftw_destroy_plan(forwardPlan);
1482 
1483  // likewise for the inverse transform
1484  fftw_plan backwardPlan = fftw_plan_dft_2d(height, width,
1485  (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(),
1486  FFTW_BACKWARD, FFTW_ESTIMATE);
1487  fftw_execute(backwardPlan);
1488  fftw_destroy_plan(backwardPlan);
1489 
1490  // do not forget to normalize the result according to the image size
1491  transformImage(srcImageRange(spatial), destImage(spatial),
1492  std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height));
1493  \endcode
1494 
1495  Note that in the creation of a plan, the height must be given
1496  first. Note also that <TT>spatial.begin()</TT> may only be passed
1497  to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the
1498  entire image. When you want to restrict operation to an ROI, you
1499  can create a copy of the ROI in an image of appropriate size, or
1500  you may use the Guru interface to FFTW.
1501 
1502  More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
1503 
1504  FFTW produces fourier images that have the DC component (the
1505  origin of the Fourier space) in the upper left corner. Often, one
1506  wants the origin in the center of the image, so that frequencies
1507  always increase towards the border of the image. This can be
1508  achieved by calling \ref moveDCToCenter(). The inverse
1509  transformation is done by \ref moveDCToUpperLeft().
1510 
1511  <b>\#include</b> <vigra/fftw3.hxx><br>
1512  Namespace: vigra
1513 */
1514 
1515 /** \addtogroup FourierTransform Fast Fourier Transform
1516 
1517  VIGRA provides a powerful C++ API for the popular <a href="http://www.fftw.org/">FFTW library</a>
1518  for fast Fourier transforms. There are two versions of the API: an older one based on image
1519  iterators (and therefore restricted to 2D) and a new one based on \ref MultiArrayView that
1520  works for arbitrary dimensions. In addition, the functions \ref convolveFFT() and
1521  \ref applyFourierFilter() provide an easy-to-use interface for FFT-based convolution,
1522  a major application of Fourier transforms.
1523 */
1524 //@{
1525 
1526 /********************************************************/
1527 /* */
1528 /* moveDCToCenter */
1529 /* */
1530 /********************************************************/
1531 
1532 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
1533  in the image center.
1534 
1535  FFTW produces fourier images where the DC component (origin of
1536  fourier space) is located in the upper left corner of the
1537  image. The quadrants are placed like this (using a 4x4 image for
1538  example):
1539 
1540  \code
1541  DC 4 3 3
1542  4 4 3 3
1543  1 1 2 2
1544  1 1 2 2
1545  \endcode
1546 
1547  After applying the function, the quadrants are at their usual places:
1548 
1549  \code
1550  2 2 1 1
1551  2 2 1 1
1552  3 3 DC 4
1553  3 3 4 4
1554  \endcode
1555 
1556  This transformation can be reversed by \ref moveDCToUpperLeft().
1557  Note that the 2D versions of this transformation must not be executed in place - input
1558  and output images must be different. In contrast, the nD version (with MultiArrayView
1559  argument) always works in-place.
1560 
1561  <b> Declarations:</b>
1562 
1563  use MultiArrayView (this works in-place, with arbitrary dimension N):
1564  \code
1565  namespace vigra {
1566  template <unsigned int N, class T, class Stride>
1567  void moveDCToCenter(MultiArrayView<N, T, Stride> a);
1568  }
1569  \endcode
1570 
1571  pass iterators explicitly (2D only, not in-place):
1572  \code
1573  namespace vigra {
1574  template <class SrcImageIterator, class SrcAccessor,
1575  class DestImageIterator, class DestAccessor>
1576  void moveDCToCenter(SrcImageIterator src_upperleft,
1577  SrcImageIterator src_lowerright, SrcAccessor sa,
1578  DestImageIterator dest_upperleft, DestAccessor da);
1579  }
1580  \endcode
1581 
1582 
1583  use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place):
1584  \code
1585  namespace vigra {
1586  template <class SrcImageIterator, class SrcAccessor,
1587  class DestImageIterator, class DestAccessor>
1588  void moveDCToCenter(
1589  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1590  pair<DestImageIterator, DestAccessor> dest);
1591  }
1592  \endcode
1593 
1594  <b> Usage:</b>
1595 
1596  <b>\#include</b> <vigra/fftw3.hxx> (for 2D variants) <br>
1597  <b>\#include</b> <vigra/multi_fft.hxx> (for nD variants) <br>
1598  Namespace: vigra
1599 
1600  \code
1601  vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
1602  ... // fill image with data
1603 
1604  // create a plan with estimated performance optimization
1605  fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
1606  (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
1607  FFTW_FORWARD, FFTW_ESTIMATE );
1608  // calculate FFT
1609  fftw_execute(forwardPlan);
1610 
1611  vigra::FFTWComplexImage rearrangedFourier(width, height);
1612  moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
1613 
1614  // delete the plan
1615  fftw_destroy_plan(forwardPlan);
1616  \endcode
1617 */
1618 doxygen_overloaded_function(template <...> void moveDCToCenter)
1619 
1620 template <class SrcImageIterator, class SrcAccessor,
1621  class DestImageIterator, class DestAccessor>
1622 void moveDCToCenter(SrcImageIterator src_upperleft,
1623  SrcImageIterator src_lowerright, SrcAccessor sa,
1624  DestImageIterator dest_upperleft, DestAccessor da)
1625 {
1626  int w = int(src_lowerright.x - src_upperleft.x);
1627  int h = int(src_lowerright.y - src_upperleft.y);
1628  int w1 = w/2;
1629  int h1 = h/2;
1630  int w2 = (w+1)/2;
1631  int h2 = (h+1)/2;
1632 
1633  // 2. Quadrant zum 4.
1634  copyImage(srcIterRange(src_upperleft,
1635  src_upperleft + Diff2D(w2, h2), sa),
1636  destIter (dest_upperleft + Diff2D(w1, h1), da));
1637 
1638  // 4. Quadrant zum 2.
1639  copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1640  src_lowerright, sa),
1641  destIter (dest_upperleft, da));
1642 
1643  // 1. Quadrant zum 3.
1644  copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1645  src_upperleft + Diff2D(w, h2), sa),
1646  destIter (dest_upperleft + Diff2D(0, h1), da));
1647 
1648  // 3. Quadrant zum 1.
1649  copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1650  src_upperleft + Diff2D(w2, h), sa),
1651  destIter (dest_upperleft + Diff2D(w1, 0), da));
1652 }
1653 
1654 template <class SrcImageIterator, class SrcAccessor,
1655  class DestImageIterator, class DestAccessor>
1656 inline void moveDCToCenter(
1657  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1658  pair<DestImageIterator, DestAccessor> dest)
1659 {
1660  moveDCToCenter(src.first, src.second, src.third,
1661  dest.first, dest.second);
1662 }
1663 
1664 /********************************************************/
1665 /* */
1666 /* moveDCToUpperLeft */
1667 /* */
1668 /********************************************************/
1669 
1670 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
1671  in the image's upper left.
1672 
1673  This function is the inversion of \ref moveDCToCenter(). See there
1674  for a detailed description and usage examples.
1675 
1676  <b> Declarations:</b>
1677 
1678  use MultiArrayView (this works in-place, with arbitrary dimension N):
1679  \code
1680  namespace vigra {
1681  template <unsigned int N, class T, class Stride>
1682  void moveDCToUpperLeft(MultiArrayView<N, T, Stride> a);
1683  }
1684  \endcode
1685 
1686  pass iterators explicitly (2D only, not in-place):
1687  \code
1688  namespace vigra {
1689  template <class SrcImageIterator, class SrcAccessor,
1690  class DestImageIterator, class DestAccessor>
1691  void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1692  SrcImageIterator src_lowerright, SrcAccessor sa,
1693  DestImageIterator dest_upperleft, DestAccessor da);
1694  }
1695  \endcode
1696 
1697 
1698  use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place):
1699  \code
1700  namespace vigra {
1701  template <class SrcImageIterator, class SrcAccessor,
1702  class DestImageIterator, class DestAccessor>
1703  void moveDCToUpperLeft(
1704  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1705  pair<DestImageIterator, DestAccessor> dest);
1706  }
1707  \endcode
1708 */
1709 doxygen_overloaded_function(template <...> void moveDCToUpperLeft)
1710 
1711 template <class SrcImageIterator, class SrcAccessor,
1712  class DestImageIterator, class DestAccessor>
1713 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1714  SrcImageIterator src_lowerright, SrcAccessor sa,
1715  DestImageIterator dest_upperleft, DestAccessor da)
1716 {
1717  int w = int(src_lowerright.x - src_upperleft.x);
1718  int h = int(src_lowerright.y - src_upperleft.y);
1719  int w2 = w/2;
1720  int h2 = h/2;
1721  int w1 = (w+1)/2;
1722  int h1 = (h+1)/2;
1723 
1724  // 2. Quadrant zum 4.
1725  copyImage(srcIterRange(src_upperleft,
1726  src_upperleft + Diff2D(w2, h2), sa),
1727  destIter (dest_upperleft + Diff2D(w1, h1), da));
1728 
1729  // 4. Quadrant zum 2.
1730  copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1731  src_lowerright, sa),
1732  destIter (dest_upperleft, da));
1733 
1734  // 1. Quadrant zum 3.
1735  copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1736  src_upperleft + Diff2D(w, h2), sa),
1737  destIter (dest_upperleft + Diff2D(0, h1), da));
1738 
1739  // 3. Quadrant zum 1.
1740  copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1741  src_upperleft + Diff2D(w2, h), sa),
1742  destIter (dest_upperleft + Diff2D(w1, 0), da));
1743 }
1744 
1745 template <class SrcImageIterator, class SrcAccessor,
1746  class DestImageIterator, class DestAccessor>
1747 inline void moveDCToUpperLeft(
1748  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1749  pair<DestImageIterator, DestAccessor> dest)
1750 {
1751  moveDCToUpperLeft(src.first, src.second, src.third,
1752  dest.first, dest.second);
1753 }
1754 
1755 namespace detail {
1756 
1757 template <class T>
1758 void
1759 fourierTransformImpl(FFTWComplexImage::const_traverser sul,
1760  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1761  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest, T sign)
1762 {
1763  int w = int(slr.x - sul.x);
1764  int h = int(slr.y - sul.y);
1765 
1766  FFTWComplexImage sworkImage, dworkImage;
1767 
1768  fftw_complex * srcPtr = (fftw_complex *)(&*sul);
1769  fftw_complex * destPtr = (fftw_complex *)(&*dul);
1770 
1771  // test for right memory layout (fftw expects a 2*width*height floats array)
1772  if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1))))
1773  {
1774  sworkImage.resize(w, h);
1775  copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
1776  srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft()));
1777  }
1778  if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1779  {
1780  dworkImage.resize(w, h);
1781  destPtr = (fftw_complex *)(&(*dworkImage.upperLeft()));
1782  }
1783 
1784  fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
1785  fftw_execute(plan);
1786  fftw_destroy_plan(plan);
1787 
1788  if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1789  {
1790  copyImage(srcImageRange(dworkImage), destIter(dul, dest));
1791  }
1792 }
1793 
1794 } // namespace detail
1795 
1796 /********************************************************/
1797 /* */
1798 /* fourierTransform */
1799 /* */
1800 /********************************************************/
1801 
1802 /** \brief Compute forward and inverse Fourier transforms.
1803 
1804  The array referring to the spatial domain (i.e. the input in a forward transform,
1805  and the output in an inverse transform) may be scalar or complex. The array representing
1806  the frequency domain (i.e. output for forward transform, input for inverse transform)
1807  must always be complex.
1808 
1809  The new implementations (those using MultiArrayView arguments) perform a normalized transform,
1810  whereas the old ones (using 2D iterators or argument objects) perform an un-normalized
1811  transform (i.e. the result of the inverse transform is scaled by the number of pixels).
1812 
1813  In general, input and output arrays must have the same shape, with the exception of the
1814  special <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C
1815  and C2R modes</a> defined by FFTW.
1816 
1817  The R2C transform reduces the redundancy in the Fourier representation of a real-valued signal:
1818  Since the Fourier representation of a real signal is symmetric, about half of the Fourier coefficients
1819  can simply be dropped. By convention, this reduction is applied to the first (innermost) dimension,
1820  such that <tt>fourier.shape(0) == spatial.shape(0)/2 + 1</tt> holds. The correct frequency domain
1821  shape can be conveniently computed by means of the function \ref fftwCorrespondingShapeR2C().
1822 
1823  Note that your program must always link against <tt>libfftw3</tt>. If you want to compute Fourier
1824  transforms for <tt>float</tt> or <tt>long double</tt> arrays, you must <i>additionally</i> link against <tt>libfftw3f</tt> and <tt>libfftw3l</tt> respectively. (Old-style functions only support <tt>double</tt>).
1825 
1826  The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a>
1827  which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e.
1828  optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning,
1829  you can use the class \ref FFTWPlan.
1830 
1831  <b> Declarations:</b>
1832 
1833  use complex-valued MultiArrayView arguments (this works for arbitrary dimension N):
1834  \code
1835  namespace vigra {
1836  template <unsigned int N, class Real, class C1, class C2>
1837  void
1838  fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1839  MultiArrayView<N, FFTWComplex<Real>, C2> out);
1840 
1841  template <unsigned int N, class Real, class C1, class C2>
1842  void
1843  fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1844  MultiArrayView<N, FFTWComplex<Real>, C2> out);
1845  }
1846  \endcode
1847 
1848  use real-valued MultiArrayView in the spatial domain, complex-valued MultiArrayView
1849  in the frequency domain (this works for arbitrary dimension N, and also supports
1850  the R2C and C2R transform, depending on the array shape in the frequency domain):
1851  \code
1852  namespace vigra {
1853  template <unsigned int N, class Real, class C1, class C2>
1854  void
1855  fourierTransform(MultiArrayView<N, Real, C1> in,
1856  MultiArrayView<N, FFTWComplex<Real>, C2> out);
1857 
1858  template <unsigned int N, class Real, class C1, class C2>
1859  void
1860  fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1861  MultiArrayView<N, Real, C2> out);
1862  }
1863  \endcode
1864 
1865  pass iterators explicitly (2D only, double only):
1866  \code
1867  namespace vigra {
1868  template <class SrcImageIterator, class SrcAccessor>
1869  void fourierTransform(SrcImageIterator srcUpperLeft,
1870  SrcImageIterator srcLowerRight, SrcAccessor src,
1871  FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest);
1872 
1873  void
1874  fourierTransformInverse(FFTWComplexImage::const_traverser sul,
1875  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1876  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
1877  }
1878  \endcode
1879 
1880  use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, double only):
1881  \code
1882  namespace vigra {
1883  template <class SrcImageIterator, class SrcAccessor>
1884  void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1885  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
1886 
1887  void
1888  fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
1889  FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
1890  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
1891  }
1892  \endcode
1893 
1894  <b> Usage:</b>
1895 
1896  <b>\#include</b> <vigra/fftw3.hxx> (old-style 2D variants)<br>
1897  <b>\#include</b> <vigra/multi_fft.hxx> (new-style nD variants)<br>
1898  Namespace: vigra
1899 
1900  old-style example (using FFTWComplexImage):
1901  \code
1902  // compute complex Fourier transform of a real image, old-style implementation
1903  vigra::DImage src(w, h);
1904  vigra::FFTWComplexImage fourier(w, h);
1905 
1906  fourierTransform(srcImageRange(src), destImage(fourier));
1907 
1908  // compute inverse Fourier transform
1909  // note that both source and destination image must be of type vigra::FFTWComplexImage
1910  vigra::FFTWComplexImage inverseFourier(w, h);
1911 
1912  fourierTransformInverse(srcImageRange(fourier), destImage(inverseFourier));
1913  \endcode
1914 
1915  new-style examples (using MultiArray):
1916  \code
1917  // compute Fourier transform of a real array, using the R2C algorithm
1918  MultiArray<2, double> src(Shape2(w, h));
1919  MultiArray<2, FFTWComplex<double> > fourier(fftwCorrespondingShapeR2C(src.shape()));
1920 
1921  fourierTransform(src, fourier);
1922 
1923  // compute inverse Fourier transform, using the C2R algorithm
1924  MultiArray<2, double> dest(src.shape());
1925  fourierTransformInverse(fourier, dest);
1926  \endcode
1927 
1928  \code
1929  // compute Fourier transform of a real array with standard algorithm
1930  MultiArray<2, double> src(Shape2(w, h));
1931  MultiArray<2, FFTWComplex<double> > fourier(src.shape());
1932 
1933  fourierTransform(src, fourier);
1934 
1935  // compute inverse Fourier transform, using the C2R algorithm
1936  MultiArray<2, double> dest(src.shape());
1937  fourierTransformInverse(fourier, dest);
1938  \endcode
1939  Complex input arrays are handled in the same way.
1940 
1941 */
1942 doxygen_overloaded_function(template <...> void fourierTransform)
1943 
1944 inline void
1945 fourierTransform(FFTWComplexImage::const_traverser sul,
1946  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1947  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
1948 {
1949  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
1950 }
1951 
1952 template <class SrcImageIterator, class SrcAccessor>
1953 void fourierTransform(SrcImageIterator srcUpperLeft,
1954  SrcImageIterator srcLowerRight, SrcAccessor sa,
1955  FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor da)
1956 {
1957  // copy real input images into a complex one...
1958  int w= srcLowerRight.x - srcUpperLeft.x;
1959  int h= srcLowerRight.y - srcUpperLeft.y;
1960 
1961  FFTWComplexImage workImage(w, h);
1962  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1963  destImage(workImage, FFTWWriteRealAccessor<>()));
1964 
1965  // ...and call the complex -> complex version of the algorithm
1966  FFTWComplexImage const & cworkImage = workImage;
1967  fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1968  destUpperLeft, da);
1969 }
1970 
1971 template <class SrcImageIterator, class SrcAccessor>
1972 inline
1973 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1974  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
1975 {
1976  fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
1977 }
1978 
1979 /** \brief Compute inverse Fourier transforms.
1980 
1981  See \ref fourierTransform() for details.
1982 */
1983 doxygen_overloaded_function(template <...> void fourierTransformInverse)
1984 
1985 inline void
1986 fourierTransformInverse(FFTWComplexImage::const_traverser sul,
1987  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1988  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
1989 {
1990  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
1991 }
1992 
1993 template <class DestImageIterator, class DestAccessor>
1994 void fourierTransformInverse(FFTWComplexImage::const_traverser sul,
1995  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1996  DestImageIterator dul, DestAccessor dest)
1997 {
1998  int w = slr.x - sul.x;
1999  int h = slr.y - sul.y;
2000 
2001  FFTWComplexImage workImage(w, h);
2002  fourierTransformInverse(sul, slr, src, workImage.upperLeft(), workImage.accessor());
2003  copyImage(srcImageRange(workImage), destIter(dul, dest));
2004 }
2005 
2006 
2007 template <class DestImageIterator, class DestAccessor>
2008 inline void
2009 fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
2010  FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
2011  pair<DestImageIterator, DestAccessor> dest)
2012 {
2013  fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second);
2014 }
2015 
2016 
2017 
2018 /********************************************************/
2019 /* */
2020 /* applyFourierFilter */
2021 /* */
2022 /********************************************************/
2023 
2024 /** \brief Apply a filter (defined in the frequency domain) to an image.
2025 
2026  After transferring the image into the frequency domain, it is
2027  multiplied pixel-wise with the filter and transformed back. The
2028  result is put into the given destination image which must have the right size.
2029  The result will be normalized to compensate for the two FFTs.
2030 
2031  If the destination image is scalar, only the real part of the result image is
2032  retained. In this case, you are responsible for choosing a filter image
2033  which ensures a zero imaginary part of the result (e.g. use a real, even symmetric
2034  filter image, or a purely imaginary, odd symmetric one).
2035 
2036  The DC entry of the filter must be in the upper left, which is the
2037  position where FFTW expects it (see \ref moveDCToUpperLeft()).
2038 
2039  See also \ref convolveFFT() for corresponding functionality on the basis of the
2040  \ref MultiArrayView interface.
2041 
2042  <b> Declarations:</b>
2043 
2044  pass 2D array views:
2045  \code
2046  namespace vigra {
2047  template <class SrcImageIterator, class SrcAccessor,
2048  class FilterImageIterator, class FilterAccessor,
2049  class DestImageIterator, class DestAccessor>
2050  void applyFourierFilter(SrcImageIterator srcUpperLeft,
2051  SrcImageIterator srcLowerRight, SrcAccessor sa,
2052  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2053  DestImageIterator destUpperLeft, DestAccessor da);
2054  }
2055  \endcode
2056 
2057  pass \ref ImageIterators and \ref DataAccessors :
2058  \code
2059  namespace vigra {
2060  template <class SrcImageIterator, class SrcAccessor,
2061  class FilterImageIterator, class FilterAccessor,
2062  class DestImageIterator, class DestAccessor>
2063  void applyFourierFilter(SrcImageIterator srcUpperLeft,
2064  SrcImageIterator srcLowerRight, SrcAccessor sa,
2065  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2066  DestImageIterator destUpperLeft, DestAccessor da);
2067  }
2068  \endcode
2069 
2070  use argument objects in conjunction with \ref ArgumentObjectFactories :
2071  \code
2072  namespace vigra {
2073  template <class SrcImageIterator, class SrcAccessor,
2074  class FilterImageIterator, class FilterAccessor,
2075  class DestImageIterator, class DestAccessor>
2076  void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2077  pair<FilterImageIterator, FilterAccessor> filter,
2078  pair<DestImageIterator, DestAccessor> dest);
2079  }
2080  \endcode
2081 
2082  <b> Usage:</b>
2083 
2084  <b>\#include</b> <vigra/fftw3.hxx><br>
2085  Namespace: vigra
2086 
2087  \code
2088  // create a Gaussian filter in Fourier space
2089  vigra::FImage gaussFilter(w, h), filter(w, h);
2090  for(int y=0; y<h; ++y)
2091  for(int x=0; x<w; ++x)
2092  {
2093  xx = float(x - w / 2) / w;
2094  yy = float(y - h / 2) / h;
2095 
2096  gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
2097  }
2098 
2099  // applyFourierFilter() expects the filter's DC in the upper left
2100  moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
2101 
2102  vigra::FFTWComplexImage result(w, h);
2103 
2104  vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
2105  \endcode
2106 
2107  For inspection of the result, \ref FFTWMagnitudeAccessor might be
2108  useful. If you want to apply the same filter repeatedly, it may be more
2109  efficient to use the FFTW functions directly with FFTW plans optimized
2110  for good performance.
2111 */
2112 doxygen_overloaded_function(template <...> void applyFourierFilter)
2113 
2114 template <class SrcImageIterator, class SrcAccessor,
2115  class FilterImageIterator, class FilterAccessor,
2116  class DestImageIterator, class DestAccessor>
2117 void applyFourierFilter(SrcImageIterator srcUpperLeft,
2118  SrcImageIterator srcLowerRight, SrcAccessor sa,
2119  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2120  DestImageIterator destUpperLeft, DestAccessor da)
2121 {
2122  // copy real input images into a complex one...
2123  int w = int(srcLowerRight.x - srcUpperLeft.x);
2124  int h = int(srcLowerRight.y - srcUpperLeft.y);
2125 
2126  FFTWComplexImage workImage(w, h);
2127  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2128  destImage(workImage, FFTWWriteRealAccessor<>()));
2129 
2130  // ...and call the impl
2131  FFTWComplexImage const & cworkImage = workImage;
2132  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2133  filterUpperLeft, fa,
2134  destUpperLeft, da);
2135 }
2136 
2137 template <class FilterImageIterator, class FilterAccessor,
2138  class DestImageIterator, class DestAccessor>
2139 inline
2140 void applyFourierFilter(
2141  FFTWComplexImage::const_traverser srcUpperLeft,
2142  FFTWComplexImage::const_traverser srcLowerRight,
2143  FFTWComplexImage::ConstAccessor sa,
2144  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2145  DestImageIterator destUpperLeft, DestAccessor da)
2146 {
2147  int w = srcLowerRight.x - srcUpperLeft.x;
2148  int h = srcLowerRight.y - srcUpperLeft.y;
2149 
2150  // test for right memory layout (fftw expects a 2*width*height floats array)
2151  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2152  applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
2153  filterUpperLeft, fa,
2154  destUpperLeft, da);
2155  else
2156  {
2157  FFTWComplexImage workImage(w, h);
2158  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2159  destImage(workImage));
2160 
2161  FFTWComplexImage const & cworkImage = workImage;
2162  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2163  filterUpperLeft, fa,
2164  destUpperLeft, da);
2165  }
2166 }
2167 
2168 template <class SrcImageIterator, class SrcAccessor,
2169  class FilterImageIterator, class FilterAccessor,
2170  class DestImageIterator, class DestAccessor>
2171 inline
2172 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2173  pair<FilterImageIterator, FilterAccessor> filter,
2174  pair<DestImageIterator, DestAccessor> dest)
2175 {
2176  applyFourierFilter(src.first, src.second, src.third,
2177  filter.first, filter.second,
2178  dest.first, dest.second);
2179 }
2180 
2181 template <class FilterImageIterator, class FilterAccessor,
2182  class DestImageIterator, class DestAccessor>
2183 void applyFourierFilterImpl(
2184  FFTWComplexImage::const_traverser srcUpperLeft,
2185  FFTWComplexImage::const_traverser srcLowerRight,
2186  FFTWComplexImage::ConstAccessor,
2187  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2188  DestImageIterator destUpperLeft, DestAccessor da)
2189 {
2190  int w = int(srcLowerRight.x - srcUpperLeft.x);
2191  int h = int(srcLowerRight.y - srcUpperLeft.y);
2192 
2193  FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
2194 
2195  // FFT from srcImage to complexResultImg
2196  fftw_plan forwardPlan=
2197  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2198  (fftw_complex *)complexResultImg.begin(),
2199  FFTW_FORWARD, FFTW_ESTIMATE );
2200  fftw_execute(forwardPlan);
2201  fftw_destroy_plan(forwardPlan);
2202 
2203  // convolve in freq. domain (in complexResultImg)
2204  combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
2205  destImage(complexResultImg), std::multiplies<FFTWComplex<> >());
2206 
2207  // FFT back into spatial domain (inplace in complexResultImg)
2208  fftw_plan backwardPlan=
2209  fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
2210  (fftw_complex *)complexResultImg.begin(),
2211  FFTW_BACKWARD, FFTW_ESTIMATE);
2212  fftw_execute(backwardPlan);
2213  fftw_destroy_plan(backwardPlan);
2214 
2215  typedef typename
2216  NumericTraits<typename DestAccessor::value_type>::isScalar
2217  isScalarResult;
2218 
2219  // normalization (after FFTs), maybe stripping imaginary part
2220  applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
2221  isScalarResult());
2222 }
2223 
2224 template <class DestImageIterator, class DestAccessor>
2225 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
2226  DestImageIterator destUpperLeft,
2227  DestAccessor da,
2228  VigraFalseType)
2229 {
2230  double normFactor= 1.0/(srcImage.width() * srcImage.height());
2231 
2232  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2233  {
2234  DestImageIterator dIt= destUpperLeft;
2235  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
2236  {
2237  da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
2238  da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
2239  }
2240  }
2241 }
2242 
2243 inline
2244 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
2245  FFTWComplexImage::traverser destUpperLeft,
2246  FFTWComplexImage::Accessor da,
2247  VigraFalseType)
2248 {
2249  transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
2250  linearIntensityTransform<FFTWComplex<> >(1.0/(srcImage.width() * srcImage.height())));
2251 }
2252 
2253 template <class DestImageIterator, class DestAccessor>
2254 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
2255  DestImageIterator destUpperLeft,
2256  DestAccessor da,
2257  VigraTrueType)
2258 {
2259  double normFactor= 1.0/(srcImage.width() * srcImage.height());
2260 
2261  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2262  {
2263  DestImageIterator dIt= destUpperLeft;
2264  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
2265  da.set(srcImage(x, y).re()*normFactor, dIt);
2266  }
2267 }
2268 
2269 /**********************************************************/
2270 /* */
2271 /* applyFourierFilterFamily */
2272 /* */
2273 /**********************************************************/
2274 
2275 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
2276 
2277  This provides the same functionality as \ref applyFourierFilter(),
2278  but applying several filters at once allows to avoid
2279  repeated Fourier transforms of the source image.
2280 
2281  Filters and result images must be stored in \ref vigra::ImageArray data
2282  structures. In contrast to \ref applyFourierFilter(), this function adjusts
2283  the size of the result images and the the length of the array.
2284 
2285  <b> Declarations:</b>
2286 
2287  pass 2D array views:
2288  \code
2289  namespace vigra {
2290  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2291  void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2292  SrcImageIterator srcLowerRight, SrcAccessor sa,
2293  const ImageArray<FilterType> &filters,
2294  ImageArray<FFTWComplexImage> &results)
2295  }
2296  \endcode
2297 
2298  pass \ref ImageIterators and \ref DataAccessors :
2299  \code
2300  namespace vigra {
2301  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2302  void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2303  SrcImageIterator srcLowerRight, SrcAccessor sa,
2304  const ImageArray<FilterType> &filters,
2305  ImageArray<FFTWComplexImage> &results)
2306  }
2307  \endcode
2308 
2309  use argument objects in conjunction with \ref ArgumentObjectFactories :
2310  \code
2311  namespace vigra {
2312  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2313  void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2314  const ImageArray<FilterType> &filters,
2315  ImageArray<FFTWComplexImage> &results)
2316  }
2317  \endcode
2318 
2319  <b> Usage:</b>
2320 
2321  <b>\#include</b> <vigra/fftw3.hxx><br>
2322  Namespace: vigra
2323 
2324  \code
2325  // assuming the presence of a real-valued image named "image" to
2326  // be filtered in this example
2327 
2328  vigra::ImageArray<vigra::FImage> filters(16, image.size());
2329 
2330  for (int i=0; i<filters.size(); i++)
2331  // create some meaningful filters here
2332  createMyFilterOfScale(i, destImage(filters[i]));
2333 
2334  vigra::ImageArray<vigra::FFTWComplexImage> results();
2335 
2336  vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
2337  \endcode
2338 */
2339 doxygen_overloaded_function(template <...> void applyFourierFilterFamily)
2340 
2341 template <class SrcImageIterator, class SrcAccessor,
2342  class FilterType, class DestImage>
2343 inline
2344 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2345  const ImageArray<FilterType> &filters,
2346  ImageArray<DestImage> &results)
2347 {
2348  applyFourierFilterFamily(src.first, src.second, src.third,
2349  filters, results);
2350 }
2351 
2352 template <class SrcImageIterator, class SrcAccessor,
2353  class FilterType, class DestImage>
2354 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2355  SrcImageIterator srcLowerRight, SrcAccessor sa,
2356  const ImageArray<FilterType> &filters,
2357  ImageArray<DestImage> &results)
2358 {
2359  int w = int(srcLowerRight.x - srcUpperLeft.x);
2360  int h = int(srcLowerRight.y - srcUpperLeft.y);
2361 
2362  FFTWComplexImage workImage(w, h);
2363  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2364  destImage(workImage, FFTWWriteRealAccessor<>()));
2365 
2366  FFTWComplexImage const & cworkImage = workImage;
2367  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2368  filters, results);
2369 }
2370 
2371 template <class FilterType, class DestImage>
2372 inline
2373 void applyFourierFilterFamily(
2374  FFTWComplexImage::const_traverser srcUpperLeft,
2375  FFTWComplexImage::const_traverser srcLowerRight,
2376  FFTWComplexImage::ConstAccessor sa,
2377  const ImageArray<FilterType> &filters,
2378  ImageArray<DestImage> &results)
2379 {
2380  int w= srcLowerRight.x - srcUpperLeft.x;
2381 
2382  // test for right memory layout (fftw expects a 2*width*height floats array)
2383  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2384  applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
2385  filters, results);
2386  else
2387  {
2388  int h = srcLowerRight.y - srcUpperLeft.y;
2389  FFTWComplexImage workImage(w, h);
2390  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2391  destImage(workImage));
2392 
2393  FFTWComplexImage const & cworkImage = workImage;
2394  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2395  filters, results);
2396  }
2397 }
2398 
2399 template <class FilterType, class DestImage>
2400 void applyFourierFilterFamilyImpl(
2401  FFTWComplexImage::const_traverser srcUpperLeft,
2402  FFTWComplexImage::const_traverser srcLowerRight,
2403  FFTWComplexImage::ConstAccessor /*sa*/,
2404  const ImageArray<FilterType> &filters,
2405  ImageArray<DestImage> &results)
2406 {
2407  // FIXME: sa is not used
2408  // (maybe check if StandardAccessor, else copy?)
2409 
2410  // make sure the filter images have the right dimensions
2411  vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
2412  "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2413 
2414  // make sure the result image array has the right dimensions
2415  results.resize(filters.size());
2416  results.resizeImages(filters.imageSize());
2417 
2418  // FFT from srcImage to freqImage
2419  int w = int(srcLowerRight.x - srcUpperLeft.x);
2420  int h = int(srcLowerRight.y - srcUpperLeft.y);
2421 
2422  FFTWComplexImage freqImage(w, h);
2423  FFTWComplexImage result(w, h);
2424 
2425  fftw_plan forwardPlan=
2426  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2427  (fftw_complex *)freqImage.begin(),
2428  FFTW_FORWARD, FFTW_ESTIMATE );
2429  fftw_execute(forwardPlan);
2430  fftw_destroy_plan(forwardPlan);
2431 
2432  fftw_plan backwardPlan=
2433  fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
2434  (fftw_complex *)result.begin(),
2435  FFTW_BACKWARD, FFTW_ESTIMATE );
2436  typedef typename
2437  NumericTraits<typename DestImage::Accessor::value_type>::isScalar
2438  isScalarResult;
2439 
2440  // convolve with filters in freq. domain
2441  for (unsigned int i= 0; i < filters.size(); i++)
2442  {
2443  combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
2444  destImage(result), std::multiplies<FFTWComplex<> >());
2445 
2446  // FFT back into spatial domain (inplace in destImage)
2447  fftw_execute(backwardPlan);
2448 
2449  // normalization (after FFTs), maybe stripping imaginary part
2450  applyFourierFilterImplNormalization(result,
2451  results[i].upperLeft(), results[i].accessor(),
2452  isScalarResult());
2453  }
2454  fftw_destroy_plan(backwardPlan);
2455 }
2456 
2457 /********************************************************/
2458 /* */
2459 /* fourierTransformReal */
2460 /* */
2461 /********************************************************/
2462 
2463 /** \brief Real Fourier transforms for even and odd boundary conditions
2464  (aka. cosine and sine transforms).
2465 
2466 
2467  If the image is real and has even symmetry, its Fourier transform
2468  is also real and has even symmetry. The Fourier transform of a real image with odd
2469  symmetry is imaginary and has odd symmetry. In either case, only about a quarter
2470  of the pixels need to be stored because the rest can be calculated from the symmetry
2471  properties. This is especially useful, if the original image is implicitly assumed
2472  to have reflective or anti-reflective boundary conditions. Then the "negative"
2473  pixel locations are defined as
2474 
2475  \code
2476  even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1)
2477  odd (anti-reflective boundary conditions): f[-1] = 0
2478  f[-x] = -f[x-2] (x = 2,...,N-1)
2479  \endcode
2480 
2481  end similar at the other boundary (see the FFTW documentation for details).
2482  This has the advantage that more efficient Fourier transforms that use only
2483  real numbers can be implemented. These are also known as cosine and sine transforms
2484  respectively.
2485 
2486  If you use the odd transform it is important to note that in the Fourier domain,
2487  the DC component is always zero and is therefore dropped from the data structure.
2488  This means that index 0 in an odd symmetric Fourier domain image refers to
2489  the <i>first</i> harmonic. This is especially important if an image is first
2490  cosine transformed (even symmetry), then in the Fourier domain multiplied
2491  with an odd symmetric filter (e.g. a first derivative) and finally transformed
2492  back to the spatial domain with a sine transform (odd symmetric). For this to work
2493  properly the image must be shifted left or up by one pixel (depending on whether
2494  the x- or y-axis is odd symmetric) before the inverse transform can be applied.
2495  (see example below).
2496 
2497  The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
2498  where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
2499  whether the x- and y-axis is to be transformed using even or odd symmetry.
2500  The same functions can be used for both the forward and inverse transforms,
2501  only the normalization changes. For signal processing, the following
2502  normalization factors are most appropriate:
2503 
2504  \code
2505  forward inverse
2506  ------------------------------------------------------------
2507  X even, Y even 1.0 4.0 * (w-1) * (h-1)
2508  X even, Y odd -1.0 -4.0 * (w-1) * (h+1)
2509  X odd, Y even -1.0 -4.0 * (w+1) * (h-1)
2510  X odd, Y odd 1.0 4.0 * (w+1) * (h+1)
2511  \endcode
2512 
2513  where <tt>w</tt> and <tt>h</tt> denote the image width and height.
2514 
2515  <b> Declarations:</b>
2516 
2517  pass 2D array views:
2518  \code
2519  namespace vigra {
2520  template <class SrcTraverser, class SrcAccessor,
2521  class DestTraverser, class DestAccessor>
2522  void
2523  fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2524  DestTraverser dul, DestAccessor dest, fftw_real norm);
2525 
2526  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2527  }
2528  \endcode
2529 
2530  pass \ref ImageIterators and \ref DataAccessors :
2531  \code
2532  namespace vigra {
2533  template <class SrcTraverser, class SrcAccessor,
2534  class DestTraverser, class DestAccessor>
2535  void
2536  fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2537  DestTraverser dul, DestAccessor dest, fftw_real norm);
2538 
2539  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2540  }
2541  \endcode
2542 
2543 
2544  use argument objects in conjunction with \ref ArgumentObjectFactories :
2545  \code
2546  namespace vigra {
2547  template <class SrcTraverser, class SrcAccessor,
2548  class DestTraverser, class DestAccessor>
2549  void
2550  fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2551  pair<DestTraverser, DestAccessor> dest, fftw_real norm);
2552 
2553  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2554  }
2555  \endcode
2556 
2557  <b> Usage:</b>
2558 
2559  <b>\#include</b> <vigra/fftw3.hxx><br>
2560  Namespace: vigra
2561 
2562  \code
2563  vigra::FImage spatial(width,height), fourier(width,height);
2564  ... // fill image with data
2565 
2566  // forward cosine transform == reflective boundary conditions
2567  fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
2568 
2569  // multiply with a first derivative of Gaussian in x-direction
2570  for(int y = 0; y < height; ++y)
2571  {
2572  for(int x = 1; x < width; ++x)
2573  {
2574  double dx = x * M_PI / (width - 1);
2575  double dy = y * M_PI / (height - 1);
2576  fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
2577  }
2578  fourier(width-1, y) = 0.0;
2579  }
2580 
2581  // inverse transform -- odd symmetry in x-direction, even in y,
2582  // due to symmetry of the filter
2583  fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
2584  (fftw_real)-4.0 * (width+1) * (height-1));
2585  \endcode
2586 */
2587 doxygen_overloaded_function(template <...> void fourierTransformReal)
2588 
2589 template <class SrcTraverser, class SrcAccessor,
2590  class DestTraverser, class DestAccessor>
2591 inline void
2592 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2593  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2594 {
2595  fourierTransformRealEE(src.first, src.second, src.third,
2596  dest.first, dest.second, norm);
2597 }
2598 
2599 template <class SrcTraverser, class SrcAccessor,
2600  class DestTraverser, class DestAccessor>
2601 inline void
2602 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2603  DestTraverser dul, DestAccessor dest, fftw_real norm)
2604 {
2605  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2606  norm, FFTW_REDFT00, FFTW_REDFT00);
2607 }
2608 
2609 template <class DestTraverser, class DestAccessor>
2610 inline void
2611 fourierTransformRealEE(
2612  FFTWRealImage::const_traverser sul,
2613  FFTWRealImage::const_traverser slr,
2614  FFTWRealImage::Accessor src,
2615  DestTraverser dul, DestAccessor dest, fftw_real norm)
2616 {
2617  int w = slr.x - sul.x;
2618 
2619  // test for right memory layout (fftw expects a width*height fftw_real array)
2620  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2621  fourierTransformRealImpl(sul, slr, dul, dest,
2622  norm, FFTW_REDFT00, FFTW_REDFT00);
2623  else
2624  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2625  norm, FFTW_REDFT00, FFTW_REDFT00);
2626 }
2627 
2628 /********************************************************************/
2629 
2630 template <class SrcTraverser, class SrcAccessor,
2631  class DestTraverser, class DestAccessor>
2632 inline void
2633 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2634  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2635 {
2636  fourierTransformRealOE(src.first, src.second, src.third,
2637  dest.first, dest.second, norm);
2638 }
2639 
2640 template <class SrcTraverser, class SrcAccessor,
2641  class DestTraverser, class DestAccessor>
2642 inline void
2643 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2644  DestTraverser dul, DestAccessor dest, fftw_real norm)
2645 {
2646  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2647  norm, FFTW_RODFT00, FFTW_REDFT00);
2648 }
2649 
2650 template <class DestTraverser, class DestAccessor>
2651 inline void
2652 fourierTransformRealOE(
2653  FFTWRealImage::const_traverser sul,
2654  FFTWRealImage::const_traverser slr,
2655  FFTWRealImage::Accessor src,
2656  DestTraverser dul, DestAccessor dest, fftw_real norm)
2657 {
2658  int w = slr.x - sul.x;
2659 
2660  // test for right memory layout (fftw expects a width*height fftw_real array)
2661  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2662  fourierTransformRealImpl(sul, slr, dul, dest,
2663  norm, FFTW_RODFT00, FFTW_REDFT00);
2664  else
2665  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2666  norm, FFTW_RODFT00, FFTW_REDFT00);
2667 }
2668 
2669 /********************************************************************/
2670 
2671 template <class SrcTraverser, class SrcAccessor,
2672  class DestTraverser, class DestAccessor>
2673 inline void
2674 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2675  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2676 {
2677  fourierTransformRealEO(src.first, src.second, src.third,
2678  dest.first, dest.second, norm);
2679 }
2680 
2681 template <class SrcTraverser, class SrcAccessor,
2682  class DestTraverser, class DestAccessor>
2683 inline void
2684 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2685  DestTraverser dul, DestAccessor dest, fftw_real norm)
2686 {
2687  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2688  norm, FFTW_REDFT00, FFTW_RODFT00);
2689 }
2690 
2691 template <class DestTraverser, class DestAccessor>
2692 inline void
2693 fourierTransformRealEO(
2694  FFTWRealImage::const_traverser sul,
2695  FFTWRealImage::const_traverser slr,
2696  FFTWRealImage::Accessor src,
2697  DestTraverser dul, DestAccessor dest, fftw_real norm)
2698 {
2699  int w = slr.x - sul.x;
2700 
2701  // test for right memory layout (fftw expects a width*height fftw_real array)
2702  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2703  fourierTransformRealImpl(sul, slr, dul, dest,
2704  norm, FFTW_REDFT00, FFTW_RODFT00);
2705  else
2706  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2707  norm, FFTW_REDFT00, FFTW_RODFT00);
2708 }
2709 
2710 /********************************************************************/
2711 
2712 template <class SrcTraverser, class SrcAccessor,
2713  class DestTraverser, class DestAccessor>
2714 inline void
2715 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2716  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2717 {
2718  fourierTransformRealOO(src.first, src.second, src.third,
2719  dest.first, dest.second, norm);
2720 }
2721 
2722 template <class SrcTraverser, class SrcAccessor,
2723  class DestTraverser, class DestAccessor>
2724 inline void
2725 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2726  DestTraverser dul, DestAccessor dest, fftw_real norm)
2727 {
2728  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2729  norm, FFTW_RODFT00, FFTW_RODFT00);
2730 }
2731 
2732 template <class DestTraverser, class DestAccessor>
2733 inline void
2734 fourierTransformRealOO(
2735  FFTWRealImage::const_traverser sul,
2736  FFTWRealImage::const_traverser slr,
2737  FFTWRealImage::Accessor src,
2738  DestTraverser dul, DestAccessor dest, fftw_real norm)
2739 {
2740  int w = slr.x - sul.x;
2741 
2742  // test for right memory layout (fftw expects a width*height fftw_real array)
2743  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2744  fourierTransformRealImpl(sul, slr, dul, dest,
2745  norm, FFTW_RODFT00, FFTW_RODFT00);
2746  else
2747  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2748  norm, FFTW_RODFT00, FFTW_RODFT00);
2749 }
2750 
2751 /*******************************************************************/
2752 
2753 template <class SrcTraverser, class SrcAccessor,
2754  class DestTraverser, class DestAccessor>
2755 void
2756 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2757  DestTraverser dul, DestAccessor dest,
2758  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2759 {
2760  FFTWRealImage workImage(slr - sul);
2761  copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2762  FFTWRealImage const & cworkImage = workImage;
2763  fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
2764  dul, dest, norm, kindx, kindy);
2765 }
2766 
2767 
2768 template <class DestTraverser, class DestAccessor>
2769 void
2770 fourierTransformRealImpl(
2771  FFTWRealImage::const_traverser sul,
2772  FFTWRealImage::const_traverser slr,
2773  DestTraverser dul, DestAccessor dest,
2774  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2775 {
2776  int w = slr.x - sul.x;
2777  int h = slr.y - sul.y;
2778  BasicImage<fftw_real> res(w, h);
2779 
2780  fftw_plan plan = fftw_plan_r2r_2d(h, w,
2781  (fftw_real *)&(*sul), (fftw_real *)res.begin(),
2782  kindy, kindx, FFTW_ESTIMATE);
2783  fftw_execute(plan);
2784  fftw_destroy_plan(plan);
2785 
2786  if(norm != 1.0)
2787  transformImage(srcImageRange(res), destIter(dul, dest),
2788  std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2789  else
2790  copyImage(srcImageRange(res), destIter(dul, dest));
2791 }
2792 
2793 
2794 //@}
2795 
2796 } // namespace vigra
2797 
2798 #endif // VIGRA_FFTW3_HXX
vigra::FFTWComplex::FFTWComplex
FFTWComplex(fftwf_complex const &o)
Definition: fftw3.hxx:202
vigra::FFTWRealAccessor
Definition: fftw3.hxx:1217
vigra::FFTWWriteRealAccessor::set
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Definition: fftw3.hxx:1339
vigra::fourierTransformInverse
void fourierTransformInverse(...)
Compute inverse Fourier transforms.
vigra::moveDCToCenter
void moveDCToCenter(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image center.
vigra::moveDCToUpperLeft
void moveDCToUpperLeft(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image's upper left...
vigra::atan2
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1654
vigra::FFTWRealImage
BasicImage< fftw_real > FFTWRealImage
Definition: fftw3.hxx:1132
vigra::FFTWComplex::value_type
Real value_type
Definition: fftw3.hxx:140
vigra::FFTWPhaseAccessor::operator()
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read phase at offset from iterator position.
Definition: fftw3.hxx:1443
vigra::IteratorTraits
Export associated information for each image iterator.
Definition: iteratortraits.hxx:109
vigra::conj
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition: fftw3.hxx:1030
vigra::imag
R imag(const FFTWComplex< R > &a)
imaginary part
Definition: fftw3.hxx:1023
vigra::FFTWComplex::FFTWComplex
FFTWComplex(value_type const &re=0.0, value_type const &im=0.0)
Definition: fftw3.hxx:169
vigra::FFTWImaginaryAccessor::set
void set(FFTWComplex< R > const &v, ITERATOR const &i) const
Write imaginary part at iterator position into a scalar.
Definition: fftw3.hxx:1300
vigra::FFTWWriteRealAccessor::set
void set(value_type const &v, ITERATOR const &i) const
Definition: fftw3.hxx:1330
vigra::FFTWComplex::SquaredNormType
value_type SquaredNormType
Definition: fftw3.hxx:164
vigra::arg
R arg(const FFTWComplex< R > &a)
phase
Definition: fftw3.hxx:1009
vigra::FFTWComplex::complex_type
FFTWReal2Complex< Real >::type complex_type
Definition: fftw3.hxx:136
vigra::applyFourierFilter
void applyFourierFilter(...)
Apply a filter (defined in the frequency domain) to an image.
vigra::FFTWSquaredMagnitudeAccessor
Definition: fftw3.hxx:1351
vigra::FFTWComplex::squaredMagnitude
SquaredNormType squaredMagnitude() const
Definition: fftw3.hxx:358
vigra::operator-
Diff2D operator-(Diff2D const &a, Diff2D const &b)
Definition: diff2d.hxx:711
vigra::FFTWRealAccessor::set
void set(FFTWComplex< R > const &v, ITERATOR const &i) const
Write real part at iterator position into a scalar.
Definition: fftw3.hxx:1250
vigra::FFTWMagnitudeAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1381
vigra::linalg::sin
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
vigra::FFTWWriteRealAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1324
vigra::FFTWWriteRealAccessor
Definition: fftw3.hxx:1319
vigra::BasicImage::Accessor
IteratorTraits< traverser >::DefaultAccessor Accessor
Definition: basicimage.hxx:573
vigra::FFTWComplex::phase
value_type phase() const
Definition: fftw3.hxx:368
vigra::FFTWRealAccessor::set
void set(FFTWComplex< R > const &v, ITERATOR const &i, DIFFERENCE d) const
Write real part at offset from iterator position into a scalar.
Definition: fftw3.hxx:1256
vigra::linalg::exp
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
vigra::applyFourierFilterFamily
void applyFourierFilterFamily(...)
Apply an array of filters (defined in the frequency domain) to an image.
vigra::FFTWComplex::operator=
FFTWComplex & operator=(fftw_complex const &o)
Definition: fftw3.hxx:255
vigra::FFTWComplex::operator=
FFTWComplex & operator=(TinyVector< T, 2 > const &o)
Definition: fftw3.hxx:310
vigra::FFTWComplex::const_iterator
value_type const * const_iterator
Definition: fftw3.hxx:156
vigra::transformImage
void transformImage(...)
Apply unary point transformation to each pixel.
vigra::FFTWSquaredMagnitudeAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1355
vigra::FFTWLogMagnitudeAccessor::operator()
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read natural log of magnitude at offset from iterator position.
Definition: fftw3.hxx:1417
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
vigra::FFTWComplex::operator=
FFTWComplex & operator=(fftwl_complex const &o)
Definition: fftw3.hxx:273
vigra::FFTWComplex::size
int size() const
Definition: fftw3.hxx:383
vigra::real
R real(const FFTWComplex< R > &a)
real part
Definition: fftw3.hxx:1016
vigra::FFTWImaginaryAccessor::set
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Write imaginary part at offset from iterator position from a scalar.
Definition: fftw3.hxx:1294
vigra::operator+
Diff2D operator+(Diff2D const &a, Diff2D const &b)
Definition: diff2d.hxx:739
vigra::FFTWImaginaryAccessor::set
void set(FFTWComplex< R > const &v, ITERATOR const &i, DIFFERENCE d) const
Write imaginary part at offset from iterator position into a scalar.
Definition: fftw3.hxx:1306
vigra::FFTWRealAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1222
vigra::FFTWComplex::operator=
FFTWComplex & operator=(std::complex< T > const &o)
Definition: fftw3.hxx:320
vigra::FFTWImaginaryAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1272
vigra::FFTWComplex::operator=
FFTWComplex & operator=(FFTWComplex const &o)
Definition: fftw3.hxx:236
vigra::FFTWSquaredMagnitudeAccessor::operator()
value_type operator()(ITERATOR const &i) const
Read squared magnitude at iterator position.
Definition: fftw3.hxx:1359
vigra::operator-=
FFTWComplex< R > & operator-=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
subtract-assignment
Definition: fftw3.hxx:867
vigra::BasicImage::traverser
BasicImageIterator< PIXELTYPE, PIXELTYPE ** > traverser
Definition: basicimage.hxx:528
vigra::FFTWSquaredMagnitudeAccessor::operator()
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read squared magnitude at offset from iterator position.
Definition: fftw3.hxx:1365
vigra::FFTWPhaseAccessor
Definition: fftw3.hxx:1429
vigra::FFTWLogMagnitudeAccessor::operator()
value_type operator()(ITERATOR const &i) const
Read natural log of magnitude at iterator position.
Definition: fftw3.hxx:1411
vigra::VectorAccessor
Accessor for items that are STL compatible vectors.
Definition: accessor.hxx:771
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
vigra::FFTWComplex::operator=
FFTWComplex & operator=(float o)
Definition: fftw3.hxx:291
vigra::BasicImageIterator
Definition: basicimage.hxx:262
vigra::FFTWRealAccessor::operator()
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read real part at offset from iterator position.
Definition: fftw3.hxx:1232
vigra::operator+=
FFTWComplex< R > & operator+=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
add-assignment
Definition: fftw3.hxx:859
vigra::linalg::log10
linalg::TemporaryMatrix< T > log10(MultiArrayView< 2, T, C > const &v)
vigra::FFTWComplex::FFTWComplex
FFTWComplex(FFTWComplex const &o)
Definition: fftw3.hxx:177
vigra::FFTWComplex::operator[]
reference operator[](int i)
Definition: fftw3.hxx:373
vigra::FFTWImaginaryAccessor::operator()
value_type operator()(ITERATOR const &i) const
Read imaginary part at iterator position.
Definition: fftw3.hxx:1276
vigra::FFTWComplex::const_reference
value_type const & const_reference
Definition: fftw3.hxx:148
vigra::operator!=
bool operator!=(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
not equal
Definition: fftw3.hxx:841
vigra::FFTWImaginaryAccessor::operator()
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read imaginary part at offset from iterator position.
Definition: fftw3.hxx:1282
vigra::FFTWComplex::FFTWComplex
FFTWComplex(FFTWComplex< U > const &o)
Definition: fftw3.hxx:186
vigra::FFTWPhaseAccessor::operator()
value_type operator()(ITERATOR const &i) const
Read phase at iterator position.
Definition: fftw3.hxx:1437
vigra::FFTWComplex::operator=
FFTWComplex & operator=(long double o)
Definition: fftw3.hxx:300
vigra::operator==
bool operator==(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
equal
Definition: fftw3.hxx:825
vigra::FFTWComplex::operator=
FFTWComplex & operator=(FFTWComplex< U > const &o)
Definition: fftw3.hxx:246
vigra::combineTwoImages
void combineTwoImages(...)
Combine two source images into destination image.
vigra::FFTWComplex::FFTWComplex
FFTWComplex(fftw_complex const &o)
Definition: fftw3.hxx:194
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::copyImage
void copyImage(...)
Copy source image into destination image.
vigra::operator*=
FFTWComplex< R > & operator*=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
multiply-assignment
Definition: fftw3.hxx:875
vigra::FFTWComplex::FFTWComplex
FFTWComplex(fftwl_complex const &o)
Definition: fftw3.hxx:210
vigra::FFTWRealAccessor::operator()
value_type operator()(ITERATOR const &i) const
Read real part at iterator position.
Definition: fftw3.hxx:1226
vigra::TinyVector
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
vigra::BasicImage::ConstAccessor
IteratorTraits< const_traverser >::DefaultAccessor ConstAccessor
Definition: basicimage.hxx:578
vigra::FFTWMagnitudeAccessor::operator()
value_type operator()(ITERATOR const &i) const
Read magnitude at iterator position.
Definition: fftw3.hxx:1385
vigra::FFTWComplex::NormType
value_type NormType
Definition: fftw3.hxx:160
vigra::FFTWImaginaryAccessor
Definition: fftw3.hxx:1268
vigra::BasicImage::const_traverser
ConstBasicImageIterator< PIXELTYPE, PIXELTYPE ** > const_traverser
Definition: basicimage.hxx:538
vigra::FFTWMagnitudeAccessor
Definition: fftw3.hxx:1377
vigra::linalg::log
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
vigra::BasicImage
Fundamental class template for images.
Definition: basicimage.hxx:475
vigra::FFTWLogMagnitudeAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1407
vigra::FFTWImaginaryAccessor::set
void set(value_type const &v, ITERATOR const &i) const
Write imaginary part at iterator position from a scalar.
Definition: fftw3.hxx:1288
vigra::FFTWComplex::operator-
FFTWComplex operator-() const
Definition: fftw3.hxx:353
vigra::fourierTransformReal
void fourierTransformReal(...)
Real Fourier transforms for even and odd boundary conditions (aka. cosine and sine transforms)...
vigra::ConstBasicImageIterator
Definition: basicimage.hxx:294
vigra::abs
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
vigra::linalg::tan
linalg::TemporaryMatrix< T > tan(MultiArrayView< 2, T, C > const &v)
vigra::FFTWComplex::iterator
value_type * iterator
Definition: fftw3.hxx:152
vigra::FFTWComplex::operator[]
const_reference operator[](int i) const
Definition: fftw3.hxx:378
vigra::operator/=
FFTWComplex< R > & operator/=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
divide-assignment
Definition: fftw3.hxx:884
vigra::FFTWComplex::magnitude
NormType magnitude() const
Definition: fftw3.hxx:363
vigra::linearIntensityTransform
LinearIntensityTransform< DestValueType, Multiplier > linearIntensityTransform(Multiplier scale, DestValueType offset)
Apply a linear transform to the source pixel values.
Definition: transformimage.hxx:800
vigra::FFTWRealAccessor::set
void set(value_type const &v, ITERATOR const &i) const
Write real part at iterator position from a scalar.
Definition: fftw3.hxx:1238
vigra::FFTWLogMagnitudeAccessor
Definition: fftw3.hxx:1403
vigra::linalg::cos
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
vigra::sign
T sign(T t)
The sign function.
Definition: mathutil.hxx:591
vigra::FFTWPhaseAccessor::value_type
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1433
vigra::FFTWMagnitudeAccessor::operator()
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read magnitude at offset from iterator position.
Definition: fftw3.hxx:1391
vigra::FFTWRealAccessor::set
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Write real part at offset from iterator position from a scalar.
Definition: fftw3.hxx:1244
vigra::FFTWComplex
Wrapper class for the FFTW complex types 'fftw_complex'.
Definition: fftw3.hxx:131
vigra::sqrt
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
vigra::FFTWComplex::reference
value_type & reference
Definition: fftw3.hxx:144
vigra::FFTWComplex::FFTWComplex
FFTWComplex(std::complex< T > const &o)
Definition: fftw3.hxx:219
vigra::FFTWComplexImage
BasicImage< FFTWComplex<> > FFTWComplexImage
Definition: fftw3.hxx:1192
vigra::FFTWComplex::FFTWComplex
FFTWComplex(TinyVector< T, 2 > const &o)
Definition: fftw3.hxx:228
vigra::FFTWComplex::operator=
FFTWComplex & operator=(double o)
Definition: fftw3.hxx:282
vigra::FFTWComplex::operator=
FFTWComplex & operator=(fftwf_complex const &o)
Definition: fftw3.hxx:264
vigra::fourierTransform
void fourierTransform(...)
Compute forward and inverse Fourier transforms.

© 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.11.1 (Fri May 19 2017)