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

fftw3.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_FFTW3_HXX
37 #define VIGRA_FFTW3_HXX
38 
39 #include <cmath>
40 #include <functional>
41 #include <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>
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  */
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  */
161 
162  /** The squared norm type (result of squaredMagnitde())
163  */
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  */
178  {
179  data_[0] = o.data_[0];
180  data_[1] = o.data_[1];
181  }
182 
183  /** Copy constructor.
184  */
185  template <class U>
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>
229  {
230  data_[0] = (Real)o[0];
231  data_[1] = (Real)o[1];
232  }
233 
234  /** Assignment.
235  */
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>
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  */
283  {
284  data_[0] = (Real)o;
285  data_[1] = 0.0;
286  return *this;
287  }
288 
289  /** Assignment.
290  */
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>
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  */
354  { return FFTWComplex(-data_[0], -data_[1]); }
355 
356  /** Squared magnitude x*conj(x)
357  */
359  { return data_[0]*data_[0]+data_[1]*data_[1]; }
360 
361  /** Magnitude (length of radius vector).
362  */
364  { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
365 
366  /** Phase angle.
367  */
369  { return VIGRA_CSTD::atan2(data_[1], data_[0]); }
370 
371  /** Access components as if number were a vector.
372  */
374  { return data_[i]; }
375 
376  /** Read components as if number were a vector.
377  */
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 
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>
860  a.re() += b.re();
861  a.im() += b.im();
862  return a;
863 }
864 
865  /// subtract-assignment
866 template <class R>
868  a.re() -= b.re();
869  a.im() -= b.im();
870  return a;
871 }
872 
873  /// multiply-assignment
874 template <class R>
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>
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>
925  a += b;
926  return a;
927 }
928 
929  /// right addition with scalar double
930 template <class R>
932  a += b;
933  return a;
934 }
935 
936  /// left addition with scalar double
937 template <class R>
939  b += a;
940  return b;
941 }
942 
943  /// subtraction
944 template <class R>
946  a -= b;
947  return a;
948 }
949 
950  /// right subtraction with scalar double
951 template <class R>
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>
1003 {
1004  return a.magnitude();
1005 }
1006 
1007  /// pahse
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>
1031 {
1032  return FFTWComplex<R>(a.re(), -a.im());
1033 }
1034 
1035  /// norm (= magnitude)
1036 template <class R>
1038 {
1039  return a.magnitude();
1040 }
1041 
1042  /// squared norm (= squared magnitude)
1043 template <class R>
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  */
1133 
1134 /********************************************************/
1135 /* */
1136 /* FFTWComplexImage */
1137 /* */
1138 /********************************************************/
1139 
1140 template<class R>
1141 struct IteratorTraits<
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  */
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>
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>
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>
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>
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>
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 phase 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>
1404 {
1405  public:
1406  /// The accessor's value type.
1407  typedef Real value_type;
1408 
1409  /// Read phase at iterator position.
1410  template <class ITERATOR>
1411  value_type operator()(ITERATOR const & i) const {
1412  return (*i).phase();
1413  }
1414 
1415  /// Read phase at offset from iterator position.
1416  template <class ITERATOR, class DIFFERENCE>
1417  value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
1418  return (i[d]).phase();
1419  }
1420 };
1421 
1422 //@}
1423 //@}
1424 
1425 /********************************************************/
1426 /* */
1427 /* Fourier Transform */
1428 /* */
1429 /********************************************************/
1430 
1431 /* \addtogroup FourierTransform Fast Fourier Transform
1432 
1433  This documentation describes the VIGRA interface to FFTW version 3. The interface
1434  to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated.
1435 
1436  VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
1437  Transform</a> package to perform Fourier transformations. VIGRA
1438  provides a wrapper for FFTW's complex number type (FFTWComplex),
1439  but FFTW's functions are used verbatim. If the image is stored as
1440  a FFTWComplexImage, the simplest call to an FFT function is like this:
1441 
1442  \code
1443  vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
1444  ... // fill image with data
1445 
1446  // create a plan with estimated performance optimization
1447  fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
1448  (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
1449  FFTW_FORWARD, FFTW_ESTIMATE );
1450  // calculate FFT (this can be repeated as often as needed,
1451  // with fresh data written into the source array)
1452  fftw_execute(forwardPlan);
1453 
1454  // release the plan memory
1455  fftw_destroy_plan(forwardPlan);
1456 
1457  // likewise for the inverse transform
1458  fftw_plan backwardPlan = fftw_plan_dft_2d(height, width,
1459  (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(),
1460  FFTW_BACKWARD, FFTW_ESTIMATE);
1461  fftw_execute(backwardPlan);
1462  fftw_destroy_plan(backwardPlan);
1463 
1464  // do not forget to normalize the result according to the image size
1465  transformImage(srcImageRange(spatial), destImage(spatial),
1466  std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height));
1467  \endcode
1468 
1469  Note that in the creation of a plan, the height must be given
1470  first. Note also that <TT>spatial.begin()</TT> may only be passed
1471  to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the
1472  entire image. When you want to restrict operation to an ROI, you
1473  can create a copy of the ROI in an image of appropriate size, or
1474  you may use the Guru interface to FFTW.
1475 
1476  More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
1477 
1478  FFTW produces fourier images that have the DC component (the
1479  origin of the Fourier space) in the upper left corner. Often, one
1480  wants the origin in the center of the image, so that frequencies
1481  always increase towards the border of the image. This can be
1482  achieved by calling \ref moveDCToCenter(). The inverse
1483  transformation is done by \ref moveDCToUpperLeft().
1484 
1485  <b>\#include</b> <vigra/fftw3.hxx><br>
1486  Namespace: vigra
1487 */
1488 
1489 /** \addtogroup FourierTransform Fast Fourier Transform
1490 
1491  VIGRA provides a powerful C++ API for the popular <a href="http://www.fftw.org/">FFTW library</a>
1492  for fast Fourier transforms. There are two versions of the API: an older one based on image
1493  iterators (and therefore restricted to 2D) and a new one based on \ref MultiArrayView that
1494  works for arbitrary dimensions. In addition, the functions \ref convolveFFT() and
1495  \ref applyFourierFilter() provide an easy-to-use interface for FFT-based convolution,
1496  a major application of Fourier transforms.
1497 */
1498 //@{
1499 
1500 /********************************************************/
1501 /* */
1502 /* moveDCToCenter */
1503 /* */
1504 /********************************************************/
1505 
1506 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
1507  in the image center.
1508 
1509  FFTW produces fourier images where the DC component (origin of
1510  fourier space) is located in the upper left corner of the
1511  image. The quadrants are placed like this (using a 4x4 image for
1512  example):
1513 
1514  \code
1515  DC 4 3 3
1516  4 4 3 3
1517  1 1 2 2
1518  1 1 2 2
1519  \endcode
1520 
1521  After applying the function, the quadrants are at their usual places:
1522 
1523  \code
1524  2 2 1 1
1525  2 2 1 1
1526  3 3 DC 4
1527  3 3 4 4
1528  \endcode
1529 
1530  This transformation can be reversed by \ref moveDCToUpperLeft().
1531  Note that the 2D versions of this transformation must not be executed in place - input
1532  and output images must be different. In contrast, the nD version (with MultiArrayView
1533  argument) always works in-place.
1534 
1535  <b> Declarations:</b>
1536 
1537  use MultiArrayView (this works in-place, with arbitrary dimension N):
1538  \code
1539  namespace vigra {
1540  template <unsigned int N, class T, class Stride>
1541  void moveDCToCenter(MultiArrayView<N, T, Stride> a);
1542  }
1543  \endcode
1544 
1545  pass iterators explicitly (2D only, not in-place):
1546  \code
1547  namespace vigra {
1548  template <class SrcImageIterator, class SrcAccessor,
1549  class DestImageIterator, class DestAccessor>
1550  void moveDCToCenter(SrcImageIterator src_upperleft,
1551  SrcImageIterator src_lowerright, SrcAccessor sa,
1552  DestImageIterator dest_upperleft, DestAccessor da);
1553  }
1554  \endcode
1555 
1556 
1557  use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place):
1558  \code
1559  namespace vigra {
1560  template <class SrcImageIterator, class SrcAccessor,
1561  class DestImageIterator, class DestAccessor>
1562  void moveDCToCenter(
1563  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1564  pair<DestImageIterator, DestAccessor> dest);
1565  }
1566  \endcode
1567 
1568  <b> Usage:</b>
1569 
1570  <b>\#include</b> <vigra/fftw3.hxx> (for 2D variants) <br>
1571  <b>\#include</b> <vigra/multi_fft.hxx> (for nD variants) <br>
1572  Namespace: vigra
1573 
1574  \code
1575  vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
1576  ... // fill image with data
1577 
1578  // create a plan with estimated performance optimization
1579  fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
1580  (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
1581  FFTW_FORWARD, FFTW_ESTIMATE );
1582  // calculate FFT
1583  fftw_execute(forwardPlan);
1584 
1585  vigra::FFTWComplexImage rearrangedFourier(width, height);
1586  moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
1587 
1588  // delete the plan
1589  fftw_destroy_plan(forwardPlan);
1590  \endcode
1591 */
1592 doxygen_overloaded_function(template <...> void moveDCToCenter)
1593 
1594 template <class SrcImageIterator, class SrcAccessor,
1595  class DestImageIterator, class DestAccessor>
1596 void moveDCToCenter(SrcImageIterator src_upperleft,
1597  SrcImageIterator src_lowerright, SrcAccessor sa,
1598  DestImageIterator dest_upperleft, DestAccessor da)
1599 {
1600  int w = int(src_lowerright.x - src_upperleft.x);
1601  int h = int(src_lowerright.y - src_upperleft.y);
1602  int w1 = w/2;
1603  int h1 = h/2;
1604  int w2 = (w+1)/2;
1605  int h2 = (h+1)/2;
1606 
1607  // 2. Quadrant zum 4.
1608  copyImage(srcIterRange(src_upperleft,
1609  src_upperleft + Diff2D(w2, h2), sa),
1610  destIter (dest_upperleft + Diff2D(w1, h1), da));
1611 
1612  // 4. Quadrant zum 2.
1613  copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1614  src_lowerright, sa),
1615  destIter (dest_upperleft, da));
1616 
1617  // 1. Quadrant zum 3.
1618  copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1619  src_upperleft + Diff2D(w, h2), sa),
1620  destIter (dest_upperleft + Diff2D(0, h1), da));
1621 
1622  // 3. Quadrant zum 1.
1623  copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1624  src_upperleft + Diff2D(w2, h), sa),
1625  destIter (dest_upperleft + Diff2D(w1, 0), da));
1626 }
1627 
1628 template <class SrcImageIterator, class SrcAccessor,
1629  class DestImageIterator, class DestAccessor>
1630 inline void moveDCToCenter(
1631  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1632  pair<DestImageIterator, DestAccessor> dest)
1633 {
1634  moveDCToCenter(src.first, src.second, src.third,
1635  dest.first, dest.second);
1636 }
1637 
1638 /********************************************************/
1639 /* */
1640 /* moveDCToUpperLeft */
1641 /* */
1642 /********************************************************/
1643 
1644 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
1645  in the image's upper left.
1646 
1647  This function is the inversion of \ref moveDCToCenter(). See there
1648  for a detailed description and usage examples.
1649 
1650  <b> Declarations:</b>
1651 
1652  use MultiArrayView (this works in-place, with arbitrary dimension N):
1653  \code
1654  namespace vigra {
1655  template <unsigned int N, class T, class Stride>
1656  void moveDCToUpperLeft(MultiArrayView<N, T, Stride> a);
1657  }
1658  \endcode
1659 
1660  pass iterators explicitly (2D only, not in-place):
1661  \code
1662  namespace vigra {
1663  template <class SrcImageIterator, class SrcAccessor,
1664  class DestImageIterator, class DestAccessor>
1665  void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1666  SrcImageIterator src_lowerright, SrcAccessor sa,
1667  DestImageIterator dest_upperleft, DestAccessor da);
1668  }
1669  \endcode
1670 
1671 
1672  use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place):
1673  \code
1674  namespace vigra {
1675  template <class SrcImageIterator, class SrcAccessor,
1676  class DestImageIterator, class DestAccessor>
1677  void moveDCToUpperLeft(
1678  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1679  pair<DestImageIterator, DestAccessor> dest);
1680  }
1681  \endcode
1682 */
1683 doxygen_overloaded_function(template <...> void moveDCToUpperLeft)
1684 
1685 template <class SrcImageIterator, class SrcAccessor,
1686  class DestImageIterator, class DestAccessor>
1687 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1688  SrcImageIterator src_lowerright, SrcAccessor sa,
1689  DestImageIterator dest_upperleft, DestAccessor da)
1690 {
1691  int w = int(src_lowerright.x - src_upperleft.x);
1692  int h = int(src_lowerright.y - src_upperleft.y);
1693  int w2 = w/2;
1694  int h2 = h/2;
1695  int w1 = (w+1)/2;
1696  int h1 = (h+1)/2;
1697 
1698  // 2. Quadrant zum 4.
1699  copyImage(srcIterRange(src_upperleft,
1700  src_upperleft + Diff2D(w2, h2), sa),
1701  destIter (dest_upperleft + Diff2D(w1, h1), da));
1702 
1703  // 4. Quadrant zum 2.
1704  copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1705  src_lowerright, sa),
1706  destIter (dest_upperleft, da));
1707 
1708  // 1. Quadrant zum 3.
1709  copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1710  src_upperleft + Diff2D(w, h2), sa),
1711  destIter (dest_upperleft + Diff2D(0, h1), da));
1712 
1713  // 3. Quadrant zum 1.
1714  copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1715  src_upperleft + Diff2D(w2, h), sa),
1716  destIter (dest_upperleft + Diff2D(w1, 0), da));
1717 }
1718 
1719 template <class SrcImageIterator, class SrcAccessor,
1720  class DestImageIterator, class DestAccessor>
1721 inline void moveDCToUpperLeft(
1722  triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1723  pair<DestImageIterator, DestAccessor> dest)
1724 {
1725  moveDCToUpperLeft(src.first, src.second, src.third,
1726  dest.first, dest.second);
1727 }
1728 
1729 namespace detail {
1730 
1731 template <class T>
1732 void
1733 fourierTransformImpl(FFTWComplexImage::const_traverser sul,
1736 {
1737  int w = int(slr.x - sul.x);
1738  int h = int(slr.y - sul.y);
1739 
1740  FFTWComplexImage sworkImage, dworkImage;
1741 
1742  fftw_complex * srcPtr = (fftw_complex *)(&*sul);
1743  fftw_complex * destPtr = (fftw_complex *)(&*dul);
1744 
1745  // test for right memory layout (fftw expects a 2*width*height floats array)
1746  if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1))))
1747  {
1748  sworkImage.resize(w, h);
1749  copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
1750  srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft()));
1751  }
1752  if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1753  {
1754  dworkImage.resize(w, h);
1755  destPtr = (fftw_complex *)(&(*dworkImage.upperLeft()));
1756  }
1757 
1758  fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
1759  fftw_execute(plan);
1760  fftw_destroy_plan(plan);
1761 
1762  if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1763  {
1764  copyImage(srcImageRange(dworkImage), destIter(dul, dest));
1765  }
1766 }
1767 
1768 } // namespace detail
1769 
1770 /********************************************************/
1771 /* */
1772 /* fourierTransform */
1773 /* */
1774 /********************************************************/
1775 
1776 /** \brief Compute forward and inverse Fourier transforms.
1777 
1778  The array referring to the spatial domain (i.e. the input in a forward transform,
1779  and the output in an inverse transform) may be scalar or complex. The array representing
1780  the frequency domain (i.e. output for forward transform, input for inverse transform)
1781  must always be complex.
1782 
1783  The new implementations (those using MultiArrayView arguments) perform a normalized transform,
1784  whereas the old ones (using 2D iterators or argument objects) perform an un-normalized
1785  transform (i.e. the result of the inverse transform is scaled by the number of pixels).
1786 
1787  In general, input and output arrays must have the same shape, with the exception of the
1788  special <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C
1789  and C2R modes</a> defined by FFTW.
1790 
1791  The R2C transform reduces the redundancy in the Fourier representation of a real-valued signal:
1792  Since the Fourier representation of a real signal is symmetric, about half of the Fourier coefficients
1793  can simply be dropped. By convention, this reduction is applied to the first (innermost) dimension,
1794  such that <tt>fourier.shape(0) == spatial.shape(0)/2 + 1</tt> holds. The correct frequency domain
1795  shape can be conveniently computed by means of the function \ref fftwCorrespondingShapeR2C().
1796 
1797  Note that your program must always link against <tt>libfftw3</tt>. If you want to compute Fourier
1798  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>).
1799 
1800  The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a>
1801  which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e.
1802  optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning,
1803  you can use the class \ref FFTWPlan.
1804 
1805  <b> Declarations:</b>
1806 
1807  use complex-valued MultiArrayView arguments (this works for arbitrary dimension N):
1808  \code
1809  namespace vigra {
1810  template <unsigned int N, class Real, class C1, class C2>
1811  void
1812  fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1813  MultiArrayView<N, FFTWComplex<Real>, C2> out);
1814 
1815  template <unsigned int N, class Real, class C1, class C2>
1816  void
1817  fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1818  MultiArrayView<N, FFTWComplex<Real>, C2> out);
1819  }
1820  \endcode
1821 
1822  use real-valued MultiArrayView in the spatial domain, complex-valued MultiArrayView
1823  in the frequency domain (this works for arbitrary dimension N, and also supports
1824  the R2C and C2R transform, depending on the array shape in the frequency domain):
1825  \code
1826  namespace vigra {
1827  template <unsigned int N, class Real, class C1, class C2>
1828  void
1829  fourierTransform(MultiArrayView<N, Real, C1> in,
1830  MultiArrayView<N, FFTWComplex<Real>, C2> out);
1831 
1832  template <unsigned int N, class Real, class C1, class C2>
1833  void
1834  fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
1835  MultiArrayView<N, Real, C2> out);
1836  }
1837  \endcode
1838 
1839  pass iterators explicitly (2D only, double only):
1840  \code
1841  namespace vigra {
1842  template <class SrcImageIterator, class SrcAccessor>
1843  void fourierTransform(SrcImageIterator srcUpperLeft,
1844  SrcImageIterator srcLowerRight, SrcAccessor src,
1845  FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest);
1846 
1847  void
1848  fourierTransformInverse(FFTWComplexImage::const_traverser sul,
1849  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
1850  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
1851  }
1852  \endcode
1853 
1854  use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, double only):
1855  \code
1856  namespace vigra {
1857  template <class SrcImageIterator, class SrcAccessor>
1858  void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1859  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
1860 
1861  void
1862  fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
1863  FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
1864  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
1865  }
1866  \endcode
1867 
1868  <b> Usage:</b>
1869 
1870  <b>\#include</b> <vigra/fftw3.hxx> (old-style 2D variants)<br>
1871  <b>\#include</b> <vigra/multi_fft.hxx> (new-style nD variants)<br>
1872  Namespace: vigra
1873 
1874  old-style example (using FFTWComplexImage):
1875  \code
1876  // compute complex Fourier transform of a real image, old-style implementation
1877  vigra::DImage src(w, h);
1878  vigra::FFTWComplexImage fourier(w, h);
1879 
1880  fourierTransform(srcImageRange(src), destImage(fourier));
1881 
1882  // compute inverse Fourier transform
1883  // note that both source and destination image must be of type vigra::FFTWComplexImage
1884  vigra::FFTWComplexImage inverseFourier(w, h);
1885 
1886  fourierTransformInverse(srcImageRange(fourier), destImage(inverseFourier));
1887  \endcode
1888 
1889  new-style examples (using MultiArray):
1890  \code
1891  // compute Fourier transform of a real array, using the R2C algorithm
1892  MultiArray<2, double> src(Shape2(w, h));
1893  MultiArray<2, FFTWComplex<double> > fourier(fftwCorrespondingShapeR2C(src.shape()));
1894 
1895  fourierTransform(src, fourier);
1896 
1897  // compute inverse Fourier transform, using the C2R algorithm
1898  MultiArray<2, double> dest(src.shape());
1899  fourierTransformInverse(fourier, dest);
1900  \endcode
1901 
1902  \code
1903  // compute Fourier transform of a real array with standard algorithm
1904  MultiArray<2, double> src(Shape2(w, h));
1905  MultiArray<2, FFTWComplex<double> > fourier(src.shape());
1906 
1907  fourierTransform(src, fourier);
1908 
1909  // compute inverse Fourier transform, using the C2R algorithm
1910  MultiArray<2, double> dest(src.shape());
1911  fourierTransformInverse(fourier, dest);
1912  \endcode
1913  Complex input arrays are handled in the same way.
1914 
1915 */
1916 doxygen_overloaded_function(template <...> void fourierTransform)
1917 
1918 inline void
1922 {
1923  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
1924 }
1925 
1926 template <class SrcImageIterator, class SrcAccessor>
1927 void fourierTransform(SrcImageIterator srcUpperLeft,
1928  SrcImageIterator srcLowerRight, SrcAccessor sa,
1930 {
1931  // copy real input images into a complex one...
1932  int w= srcLowerRight.x - srcUpperLeft.x;
1933  int h= srcLowerRight.y - srcUpperLeft.y;
1934 
1935  FFTWComplexImage workImage(w, h);
1936  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1937  destImage(workImage, FFTWWriteRealAccessor<>()));
1938 
1939  // ...and call the complex -> complex version of the algorithm
1940  FFTWComplexImage const & cworkImage = workImage;
1941  fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1942  destUpperLeft, da);
1943 }
1944 
1945 template <class SrcImageIterator, class SrcAccessor>
1946 inline
1947 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1948  pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
1949 {
1950  fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
1951 }
1952 
1953 /** \brief Compute inverse Fourier transforms.
1954 
1955  See \ref fourierTransform() for details.
1956 */
1957 doxygen_overloaded_function(template <...> void fourierTransformInverse)
1958 
1959 inline void
1963 {
1964  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
1965 }
1966 
1967 template <class DestImageIterator, class DestAccessor>
1970  DestImageIterator dul, DestAccessor dest)
1971 {
1972  int w = slr.x - sul.x;
1973  int h = slr.y - sul.y;
1974 
1975  FFTWComplexImage workImage(w, h);
1976  fourierTransformInverse(sul, slr, src, workImage.upperLeft(), workImage.accessor());
1977  copyImage(srcImageRange(workImage), destIter(dul, dest));
1978 }
1979 
1980 
1981 template <class DestImageIterator, class DestAccessor>
1982 inline void
1985  pair<DestImageIterator, DestAccessor> dest)
1986 {
1987  fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second);
1988 }
1989 
1990 
1991 
1992 /********************************************************/
1993 /* */
1994 /* applyFourierFilter */
1995 /* */
1996 /********************************************************/
1997 
1998 /** \brief Apply a filter (defined in the frequency domain) to an image.
1999 
2000  After transferring the image into the frequency domain, it is
2001  multiplied pixel-wise with the filter and transformed back. The
2002  result is put into the given destination image which must have the right size.
2003  The result will be normalized to compensate for the two FFTs.
2004 
2005  If the destination image is scalar, only the real part of the result image is
2006  retained. In this case, you are responsible for choosing a filter image
2007  which ensures a zero imaginary part of the result (e.g. use a real, even symmetric
2008  filter image, or a purely imaginary, odd symmetric one).
2009 
2010  The DC entry of the filter must be in the upper left, which is the
2011  position where FFTW expects it (see \ref moveDCToUpperLeft()).
2012 
2013  See also \ref convolveFFT() for corresponding functionality on the basis of the
2014  \ref MultiArrayView interface.
2015 
2016  <b> Declarations:</b>
2017 
2018  pass arguments explicitly:
2019  \code
2020  namespace vigra {
2021  template <class SrcImageIterator, class SrcAccessor,
2022  class FilterImageIterator, class FilterAccessor,
2023  class DestImageIterator, class DestAccessor>
2024  void applyFourierFilter(SrcImageIterator srcUpperLeft,
2025  SrcImageIterator srcLowerRight, SrcAccessor sa,
2026  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2027  DestImageIterator destUpperLeft, DestAccessor da);
2028  }
2029  \endcode
2030 
2031  use argument objects in conjunction with \ref ArgumentObjectFactories :
2032  \code
2033  namespace vigra {
2034  template <class SrcImageIterator, class SrcAccessor,
2035  class FilterImageIterator, class FilterAccessor,
2036  class DestImageIterator, class DestAccessor>
2037  void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2038  pair<FilterImageIterator, FilterAccessor> filter,
2039  pair<DestImageIterator, DestAccessor> dest);
2040  }
2041  \endcode
2042 
2043  <b> Usage:</b>
2044 
2045  <b>\#include</b> <vigra/fftw3.hxx><br>
2046  Namespace: vigra
2047 
2048  \code
2049  // create a Gaussian filter in Fourier space
2050  vigra::FImage gaussFilter(w, h), filter(w, h);
2051  for(int y=0; y<h; ++y)
2052  for(int x=0; x<w; ++x)
2053  {
2054  xx = float(x - w / 2) / w;
2055  yy = float(y - h / 2) / h;
2056 
2057  gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
2058  }
2059 
2060  // applyFourierFilter() expects the filter's DC in the upper left
2061  moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
2062 
2063  vigra::FFTWComplexImage result(w, h);
2064 
2065  vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
2066  \endcode
2067 
2068  For inspection of the result, \ref FFTWMagnitudeAccessor might be
2069  useful. If you want to apply the same filter repeatedly, it may be more
2070  efficient to use the FFTW functions directly with FFTW plans optimized
2071  for good performance.
2072 */
2073 doxygen_overloaded_function(template <...> void applyFourierFilter)
2074 
2075 template <class SrcImageIterator, class SrcAccessor,
2076  class FilterImageIterator, class FilterAccessor,
2077  class DestImageIterator, class DestAccessor>
2078 void applyFourierFilter(SrcImageIterator srcUpperLeft,
2079  SrcImageIterator srcLowerRight, SrcAccessor sa,
2080  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2081  DestImageIterator destUpperLeft, DestAccessor da)
2082 {
2083  // copy real input images into a complex one...
2084  int w = int(srcLowerRight.x - srcUpperLeft.x);
2085  int h = int(srcLowerRight.y - srcUpperLeft.y);
2086 
2087  FFTWComplexImage workImage(w, h);
2088  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2089  destImage(workImage, FFTWWriteRealAccessor<>()));
2090 
2091  // ...and call the impl
2092  FFTWComplexImage const & cworkImage = workImage;
2093  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2094  filterUpperLeft, fa,
2095  destUpperLeft, da);
2096 }
2097 
2098 template <class FilterImageIterator, class FilterAccessor,
2099  class DestImageIterator, class DestAccessor>
2100 inline
2101 void applyFourierFilter(
2102  FFTWComplexImage::const_traverser srcUpperLeft,
2103  FFTWComplexImage::const_traverser srcLowerRight,
2105  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2106  DestImageIterator destUpperLeft, DestAccessor da)
2107 {
2108  int w = srcLowerRight.x - srcUpperLeft.x;
2109  int h = srcLowerRight.y - srcUpperLeft.y;
2110 
2111  // test for right memory layout (fftw expects a 2*width*height floats array)
2112  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2113  applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
2114  filterUpperLeft, fa,
2115  destUpperLeft, da);
2116  else
2117  {
2118  FFTWComplexImage workImage(w, h);
2119  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2120  destImage(workImage));
2121 
2122  FFTWComplexImage const & cworkImage = workImage;
2123  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2124  filterUpperLeft, fa,
2125  destUpperLeft, da);
2126  }
2127 }
2128 
2129 template <class SrcImageIterator, class SrcAccessor,
2130  class FilterImageIterator, class FilterAccessor,
2131  class DestImageIterator, class DestAccessor>
2132 inline
2133 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2134  pair<FilterImageIterator, FilterAccessor> filter,
2135  pair<DestImageIterator, DestAccessor> dest)
2136 {
2137  applyFourierFilter(src.first, src.second, src.third,
2138  filter.first, filter.second,
2139  dest.first, dest.second);
2140 }
2141 
2142 template <class FilterImageIterator, class FilterAccessor,
2143  class DestImageIterator, class DestAccessor>
2144 void applyFourierFilterImpl(
2145  FFTWComplexImage::const_traverser srcUpperLeft,
2146  FFTWComplexImage::const_traverser srcLowerRight,
2148  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2149  DestImageIterator destUpperLeft, DestAccessor da)
2150 {
2151  int w = int(srcLowerRight.x - srcUpperLeft.x);
2152  int h = int(srcLowerRight.y - srcUpperLeft.y);
2153 
2154  FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
2155 
2156  // FFT from srcImage to complexResultImg
2157  fftw_plan forwardPlan=
2158  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2159  (fftw_complex *)complexResultImg.begin(),
2160  FFTW_FORWARD, FFTW_ESTIMATE );
2161  fftw_execute(forwardPlan);
2162  fftw_destroy_plan(forwardPlan);
2163 
2164  // convolve in freq. domain (in complexResultImg)
2165  combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
2166  destImage(complexResultImg), std::multiplies<FFTWComplex<> >());
2167 
2168  // FFT back into spatial domain (inplace in complexResultImg)
2169  fftw_plan backwardPlan=
2170  fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
2171  (fftw_complex *)complexResultImg.begin(),
2172  FFTW_BACKWARD, FFTW_ESTIMATE);
2173  fftw_execute(backwardPlan);
2174  fftw_destroy_plan(backwardPlan);
2175 
2176  typedef typename
2177  NumericTraits<typename DestAccessor::value_type>::isScalar
2178  isScalarResult;
2179 
2180  // normalization (after FFTs), maybe stripping imaginary part
2181  applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
2182  isScalarResult());
2183 }
2184 
2185 template <class DestImageIterator, class DestAccessor>
2186 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
2187  DestImageIterator destUpperLeft,
2188  DestAccessor da,
2189  VigraFalseType)
2190 {
2191  double normFactor= 1.0/(srcImage.width() * srcImage.height());
2192 
2193  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2194  {
2195  DestImageIterator dIt= destUpperLeft;
2196  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
2197  {
2198  da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
2199  da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
2200  }
2201  }
2202 }
2203 
2204 inline
2205 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
2206  FFTWComplexImage::traverser destUpperLeft,
2208  VigraFalseType)
2209 {
2210  transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
2211  linearIntensityTransform<FFTWComplex<> >(1.0/(srcImage.width() * srcImage.height())));
2212 }
2213 
2214 template <class DestImageIterator, class DestAccessor>
2215 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
2216  DestImageIterator destUpperLeft,
2217  DestAccessor da,
2218  VigraTrueType)
2219 {
2220  double normFactor= 1.0/(srcImage.width() * srcImage.height());
2221 
2222  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2223  {
2224  DestImageIterator dIt= destUpperLeft;
2225  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
2226  da.set(srcImage(x, y).re()*normFactor, dIt);
2227  }
2228 }
2229 
2230 /**********************************************************/
2231 /* */
2232 /* applyFourierFilterFamily */
2233 /* */
2234 /**********************************************************/
2235 
2236 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
2237 
2238  This provides the same functionality as \ref applyFourierFilter(),
2239  but applying several filters at once allows to avoid
2240  repeated Fourier transforms of the source image.
2241 
2242  Filters and result images must be stored in \ref vigra::ImageArray data
2243  structures. In contrast to \ref applyFourierFilter(), this function adjusts
2244  the size of the result images and the the length of the array.
2245 
2246  <b> Declarations:</b>
2247 
2248  pass arguments explicitly:
2249  \code
2250  namespace vigra {
2251  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2252  void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2253  SrcImageIterator srcLowerRight, SrcAccessor sa,
2254  const ImageArray<FilterType> &filters,
2255  ImageArray<FFTWComplexImage> &results)
2256  }
2257  \endcode
2258 
2259  use argument objects in conjunction with \ref ArgumentObjectFactories :
2260  \code
2261  namespace vigra {
2262  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2263  void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2264  const ImageArray<FilterType> &filters,
2265  ImageArray<FFTWComplexImage> &results)
2266  }
2267  \endcode
2268 
2269  <b> Usage:</b>
2270 
2271  <b>\#include</b> <vigra/fftw3.hxx><br>
2272  Namespace: vigra
2273 
2274  \code
2275  // assuming the presence of a real-valued image named "image" to
2276  // be filtered in this example
2277 
2278  vigra::ImageArray<vigra::FImage> filters(16, image.size());
2279 
2280  for (int i=0; i<filters.size(); i++)
2281  // create some meaningful filters here
2282  createMyFilterOfScale(i, destImage(filters[i]));
2283 
2284  vigra::ImageArray<vigra::FFTWComplexImage> results();
2285 
2286  vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
2287  \endcode
2288 */
2289 doxygen_overloaded_function(template <...> void applyFourierFilterFamily)
2290 
2291 template <class SrcImageIterator, class SrcAccessor,
2292  class FilterType, class DestImage>
2293 inline
2294 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2295  const ImageArray<FilterType> &filters,
2296  ImageArray<DestImage> &results)
2297 {
2298  applyFourierFilterFamily(src.first, src.second, src.third,
2299  filters, results);
2300 }
2301 
2302 template <class SrcImageIterator, class SrcAccessor,
2303  class FilterType, class DestImage>
2304 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2305  SrcImageIterator srcLowerRight, SrcAccessor sa,
2306  const ImageArray<FilterType> &filters,
2307  ImageArray<DestImage> &results)
2308 {
2309  int w = int(srcLowerRight.x - srcUpperLeft.x);
2310  int h = int(srcLowerRight.y - srcUpperLeft.y);
2311 
2312  FFTWComplexImage workImage(w, h);
2313  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2314  destImage(workImage, FFTWWriteRealAccessor<>()));
2315 
2316  FFTWComplexImage const & cworkImage = workImage;
2317  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2318  filters, results);
2319 }
2320 
2321 template <class FilterType, class DestImage>
2322 inline
2324  FFTWComplexImage::const_traverser srcUpperLeft,
2325  FFTWComplexImage::const_traverser srcLowerRight,
2327  const ImageArray<FilterType> &filters,
2328  ImageArray<DestImage> &results)
2329 {
2330  int w= srcLowerRight.x - srcUpperLeft.x;
2331 
2332  // test for right memory layout (fftw expects a 2*width*height floats array)
2333  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2334  applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
2335  filters, results);
2336  else
2337  {
2338  int h = srcLowerRight.y - srcUpperLeft.y;
2339  FFTWComplexImage workImage(w, h);
2340  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2341  destImage(workImage));
2342 
2343  FFTWComplexImage const & cworkImage = workImage;
2344  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2345  filters, results);
2346  }
2347 }
2348 
2349 template <class FilterType, class DestImage>
2350 void applyFourierFilterFamilyImpl(
2351  FFTWComplexImage::const_traverser srcUpperLeft,
2352  FFTWComplexImage::const_traverser srcLowerRight,
2354  const ImageArray<FilterType> &filters,
2355  ImageArray<DestImage> &results)
2356 {
2357  // FIXME: sa is not used
2358  // (maybe check if StandardAccessor, else copy?)
2359 
2360  // make sure the filter images have the right dimensions
2361  vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
2362  "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2363 
2364  // make sure the result image array has the right dimensions
2365  results.resize(filters.size());
2366  results.resizeImages(filters.imageSize());
2367 
2368  // FFT from srcImage to freqImage
2369  int w = int(srcLowerRight.x - srcUpperLeft.x);
2370  int h = int(srcLowerRight.y - srcUpperLeft.y);
2371 
2372  FFTWComplexImage freqImage(w, h);
2373  FFTWComplexImage result(w, h);
2374 
2375  fftw_plan forwardPlan=
2376  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2377  (fftw_complex *)freqImage.begin(),
2378  FFTW_FORWARD, FFTW_ESTIMATE );
2379  fftw_execute(forwardPlan);
2380  fftw_destroy_plan(forwardPlan);
2381 
2382  fftw_plan backwardPlan=
2383  fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
2384  (fftw_complex *)result.begin(),
2385  FFTW_BACKWARD, FFTW_ESTIMATE );
2386  typedef typename
2387  NumericTraits<typename DestImage::Accessor::value_type>::isScalar
2388  isScalarResult;
2389 
2390  // convolve with filters in freq. domain
2391  for (unsigned int i= 0; i < filters.size(); i++)
2392  {
2393  combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
2394  destImage(result), std::multiplies<FFTWComplex<> >());
2395 
2396  // FFT back into spatial domain (inplace in destImage)
2397  fftw_execute(backwardPlan);
2398 
2399  // normalization (after FFTs), maybe stripping imaginary part
2400  applyFourierFilterImplNormalization(result,
2401  results[i].upperLeft(), results[i].accessor(),
2402  isScalarResult());
2403  }
2404  fftw_destroy_plan(backwardPlan);
2405 }
2406 
2407 /********************************************************/
2408 /* */
2409 /* fourierTransformReal */
2410 /* */
2411 /********************************************************/
2412 
2413 /** \brief Real Fourier transforms for even and odd boundary conditions
2414  (aka. cosine and sine transforms).
2415 
2416 
2417  If the image is real and has even symmetry, its Fourier transform
2418  is also real and has even symmetry. The Fourier transform of a real image with odd
2419  symmetry is imaginary and has odd symmetry. In either case, only about a quarter
2420  of the pixels need to be stored because the rest can be calculated from the symmetry
2421  properties. This is especially useful, if the original image is implicitly assumed
2422  to have reflective or anti-reflective boundary conditions. Then the "negative"
2423  pixel locations are defined as
2424 
2425  \code
2426  even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1)
2427  odd (anti-reflective boundary conditions): f[-1] = 0
2428  f[-x] = -f[x-2] (x = 2,...,N-1)
2429  \endcode
2430 
2431  end similar at the other boundary (see the FFTW documentation for details).
2432  This has the advantage that more efficient Fourier transforms that use only
2433  real numbers can be implemented. These are also known as cosine and sine transforms
2434  respectively.
2435 
2436  If you use the odd transform it is important to note that in the Fourier domain,
2437  the DC component is always zero and is therefore dropped from the data structure.
2438  This means that index 0 in an odd symmetric Fourier domain image refers to
2439  the <i>first</i> harmonic. This is especially important if an image is first
2440  cosine transformed (even symmetry), then in the Fourier domain multiplied
2441  with an odd symmetric filter (e.g. a first derivative) and finally transformed
2442  back to the spatial domain with a sine transform (odd symmetric). For this to work
2443  properly the image must be shifted left or up by one pixel (depending on whether
2444  the x- or y-axis is odd symmetric) before the inverse transform can be applied.
2445  (see example below).
2446 
2447  The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
2448  where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
2449  whether the x- and y-axis is to be transformed using even or odd symmetry.
2450  The same functions can be used for both the forward and inverse transforms,
2451  only the normalization changes. For signal processing, the following
2452  normalization factors are most appropriate:
2453 
2454  \code
2455  forward inverse
2456  ------------------------------------------------------------
2457  X even, Y even 1.0 4.0 * (w-1) * (h-1)
2458  X even, Y odd -1.0 -4.0 * (w-1) * (h+1)
2459  X odd, Y even -1.0 -4.0 * (w+1) * (h-1)
2460  X odd, Y odd 1.0 4.0 * (w+1) * (h+1)
2461  \endcode
2462 
2463  where <tt>w</tt> and <tt>h</tt> denote the image width and height.
2464 
2465  <b> Declarations:</b>
2466 
2467  pass arguments explicitly:
2468  \code
2469  namespace vigra {
2470  template <class SrcTraverser, class SrcAccessor,
2471  class DestTraverser, class DestAccessor>
2472  void
2473  fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2474  DestTraverser dul, DestAccessor dest, fftw_real norm);
2475 
2476  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2477  }
2478  \endcode
2479 
2480 
2481  use argument objects in conjunction with \ref ArgumentObjectFactories :
2482  \code
2483  namespace vigra {
2484  template <class SrcTraverser, class SrcAccessor,
2485  class DestTraverser, class DestAccessor>
2486  void
2487  fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2488  pair<DestTraverser, DestAccessor> dest, fftw_real norm);
2489 
2490  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2491  }
2492  \endcode
2493 
2494  <b> Usage:</b>
2495 
2496  <b>\#include</b> <vigra/fftw3.hxx><br>
2497  Namespace: vigra
2498 
2499  \code
2500  vigra::FImage spatial(width,height), fourier(width,height);
2501  ... // fill image with data
2502 
2503  // forward cosine transform == reflective boundary conditions
2504  fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
2505 
2506  // multiply with a first derivative of Gaussian in x-direction
2507  for(int y = 0; y < height; ++y)
2508  {
2509  for(int x = 1; x < width; ++x)
2510  {
2511  double dx = x * M_PI / (width - 1);
2512  double dy = y * M_PI / (height - 1);
2513  fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
2514  }
2515  fourier(width-1, y) = 0.0;
2516  }
2517 
2518  // inverse transform -- odd symmetry in x-direction, even in y,
2519  // due to symmetry of the filter
2520  fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
2521  (fftw_real)-4.0 * (width+1) * (height-1));
2522  \endcode
2523 */
2524 doxygen_overloaded_function(template <...> void fourierTransformReal)
2525 
2526 template <class SrcTraverser, class SrcAccessor,
2527  class DestTraverser, class DestAccessor>
2528 inline void
2529 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2530  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2531 {
2532  fourierTransformRealEE(src.first, src.second, src.third,
2533  dest.first, dest.second, norm);
2534 }
2535 
2536 template <class SrcTraverser, class SrcAccessor,
2537  class DestTraverser, class DestAccessor>
2538 inline void
2539 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2540  DestTraverser dul, DestAccessor dest, fftw_real norm)
2541 {
2542  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2543  norm, FFTW_REDFT00, FFTW_REDFT00);
2544 }
2545 
2546 template <class DestTraverser, class DestAccessor>
2547 inline void
2548 fourierTransformRealEE(
2552  DestTraverser dul, DestAccessor dest, fftw_real norm)
2553 {
2554  int w = slr.x - sul.x;
2555 
2556  // test for right memory layout (fftw expects a width*height fftw_real array)
2557  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2558  fourierTransformRealImpl(sul, slr, dul, dest,
2559  norm, FFTW_REDFT00, FFTW_REDFT00);
2560  else
2561  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2562  norm, FFTW_REDFT00, FFTW_REDFT00);
2563 }
2564 
2565 /********************************************************************/
2566 
2567 template <class SrcTraverser, class SrcAccessor,
2568  class DestTraverser, class DestAccessor>
2569 inline void
2570 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2571  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2572 {
2573  fourierTransformRealOE(src.first, src.second, src.third,
2574  dest.first, dest.second, norm);
2575 }
2576 
2577 template <class SrcTraverser, class SrcAccessor,
2578  class DestTraverser, class DestAccessor>
2579 inline void
2580 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2581  DestTraverser dul, DestAccessor dest, fftw_real norm)
2582 {
2583  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2584  norm, FFTW_RODFT00, FFTW_REDFT00);
2585 }
2586 
2587 template <class DestTraverser, class DestAccessor>
2588 inline void
2589 fourierTransformRealOE(
2593  DestTraverser dul, DestAccessor dest, fftw_real norm)
2594 {
2595  int w = slr.x - sul.x;
2596 
2597  // test for right memory layout (fftw expects a width*height fftw_real array)
2598  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2599  fourierTransformRealImpl(sul, slr, dul, dest,
2600  norm, FFTW_RODFT00, FFTW_REDFT00);
2601  else
2602  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2603  norm, FFTW_RODFT00, FFTW_REDFT00);
2604 }
2605 
2606 /********************************************************************/
2607 
2608 template <class SrcTraverser, class SrcAccessor,
2609  class DestTraverser, class DestAccessor>
2610 inline void
2611 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2612  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2613 {
2614  fourierTransformRealEO(src.first, src.second, src.third,
2615  dest.first, dest.second, norm);
2616 }
2617 
2618 template <class SrcTraverser, class SrcAccessor,
2619  class DestTraverser, class DestAccessor>
2620 inline void
2621 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2622  DestTraverser dul, DestAccessor dest, fftw_real norm)
2623 {
2624  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2625  norm, FFTW_REDFT00, FFTW_RODFT00);
2626 }
2627 
2628 template <class DestTraverser, class DestAccessor>
2629 inline void
2630 fourierTransformRealEO(
2634  DestTraverser dul, DestAccessor dest, fftw_real norm)
2635 {
2636  int w = slr.x - sul.x;
2637 
2638  // test for right memory layout (fftw expects a width*height fftw_real array)
2639  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2640  fourierTransformRealImpl(sul, slr, dul, dest,
2641  norm, FFTW_REDFT00, FFTW_RODFT00);
2642  else
2643  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2644  norm, FFTW_REDFT00, FFTW_RODFT00);
2645 }
2646 
2647 /********************************************************************/
2648 
2649 template <class SrcTraverser, class SrcAccessor,
2650  class DestTraverser, class DestAccessor>
2651 inline void
2652 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2653  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2654 {
2655  fourierTransformRealOO(src.first, src.second, src.third,
2656  dest.first, dest.second, norm);
2657 }
2658 
2659 template <class SrcTraverser, class SrcAccessor,
2660  class DestTraverser, class DestAccessor>
2661 inline void
2662 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2663  DestTraverser dul, DestAccessor dest, fftw_real norm)
2664 {
2665  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2666  norm, FFTW_RODFT00, FFTW_RODFT00);
2667 }
2668 
2669 template <class DestTraverser, class DestAccessor>
2670 inline void
2671 fourierTransformRealOO(
2675  DestTraverser dul, DestAccessor dest, fftw_real norm)
2676 {
2677  int w = slr.x - sul.x;
2678 
2679  // test for right memory layout (fftw expects a width*height fftw_real array)
2680  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2681  fourierTransformRealImpl(sul, slr, dul, dest,
2682  norm, FFTW_RODFT00, FFTW_RODFT00);
2683  else
2684  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2685  norm, FFTW_RODFT00, FFTW_RODFT00);
2686 }
2687 
2688 /*******************************************************************/
2689 
2690 template <class SrcTraverser, class SrcAccessor,
2691  class DestTraverser, class DestAccessor>
2692 void
2693 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2694  DestTraverser dul, DestAccessor dest,
2695  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2696 {
2697  FFTWRealImage workImage(slr - sul);
2698  copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2699  FFTWRealImage const & cworkImage = workImage;
2700  fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
2701  dul, dest, norm, kindx, kindy);
2702 }
2703 
2704 
2705 template <class DestTraverser, class DestAccessor>
2706 void
2707 fourierTransformRealImpl(
2710  DestTraverser dul, DestAccessor dest,
2711  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2712 {
2713  int w = slr.x - sul.x;
2714  int h = slr.y - sul.y;
2715  BasicImage<fftw_real> res(w, h);
2716 
2717  fftw_plan plan = fftw_plan_r2r_2d(h, w,
2718  (fftw_real *)&(*sul), (fftw_real *)res.begin(),
2719  kindy, kindx, FFTW_ESTIMATE);
2720  fftw_execute(plan);
2721  fftw_destroy_plan(plan);
2722 
2723  if(norm != 1.0)
2724  transformImage(srcImageRange(res), destIter(dul, dest),
2725  std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2726  else
2727  copyImage(srcImageRange(res), destIter(dul, dest));
2728 }
2729 
2730 
2731 //@}
2732 
2733 } // namespace vigra
2734 
2735 #endif // VIGRA_FFTW3_HXX
FFTWComplex(fftwf_complex const &o)
Definition: fftw3.hxx:202
Definition: fftw3.hxx:1217
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Definition: fftw3.hxx:1339
void fourierTransformInverse(...)
Compute inverse Fourier transforms.
void moveDCToCenter(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image center.
void moveDCToUpperLeft(...)
Rearrange the quadrants of a Fourier image so that the origin is in the image's upper left...
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:1640
BasicImage< fftw_real > FFTWRealImage
Definition: fftw3.hxx:1132
Real value_type
Definition: fftw3.hxx:140
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read phase at offset from iterator position.
Definition: fftw3.hxx:1417
Export associated information for each image iterator.
Definition: iteratortraits.hxx:109
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition: fftw3.hxx:1030
R imag(const FFTWComplex< R > &a)
imaginary part
Definition: fftw3.hxx:1023
FFTWComplex(value_type const &re=0.0, value_type const &im=0.0)
Definition: fftw3.hxx:169
void set(FFTWComplex< R > const &v, ITERATOR const &i) const
Write imaginary part at iterator position into a scalar.
Definition: fftw3.hxx:1300
void set(value_type const &v, ITERATOR const &i) const
Definition: fftw3.hxx:1330
value_type SquaredNormType
Definition: fftw3.hxx:164
R arg(const FFTWComplex< R > &a)
pahse
Definition: fftw3.hxx:1009
FFTWReal2Complex< Real >::type complex_type
Definition: fftw3.hxx:136
void applyFourierFilter(...)
Apply a filter (defined in the frequency domain) to an image.
Definition: fftw3.hxx:1351
SquaredNormType squaredMagnitude() const
Definition: fftw3.hxx:358
Diff2D operator-(Diff2D const &a, Diff2D const &b)
Definition: diff2d.hxx:711
void set(FFTWComplex< R > const &v, ITERATOR const &i) const
Write real part at iterator position into a scalar.
Definition: fftw3.hxx:1250
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1381
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1324
Definition: fftw3.hxx:1319
IteratorTraits< traverser >::DefaultAccessor Accessor
Definition: basicimage.hxx:571
value_type phase() const
Definition: fftw3.hxx:368
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
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void applyFourierFilterFamily(...)
Apply an array of filters (defined in the frequency domain) to an image.
FFTWComplex & operator=(fftw_complex const &o)
Definition: fftw3.hxx:255
FFTWComplex & operator=(TinyVector< T, 2 > const &o)
Definition: fftw3.hxx:310
value_type const * const_iterator
Definition: fftw3.hxx:156
void transformImage(...)
Apply unary point transformation to each pixel.
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1355
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
FFTWComplex & operator=(fftwl_complex const &o)
Definition: fftw3.hxx:273
int size() const
Definition: fftw3.hxx:383
R real(const FFTWComplex< R > &a)
real part
Definition: fftw3.hxx:1016
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
Diff2D operator+(Diff2D const &a, Diff2D const &b)
Definition: diff2d.hxx:739
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
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1222
FFTWComplex & operator=(std::complex< T > const &o)
Definition: fftw3.hxx:320
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1272
FFTWComplex & operator=(FFTWComplex const &o)
Definition: fftw3.hxx:236
value_type operator()(ITERATOR const &i) const
Read squared magnitude at iterator position.
Definition: fftw3.hxx:1359
FFTWComplex< R > & operator-=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
subtract-assignment
Definition: fftw3.hxx:867
BasicImageIterator< PIXELTYPE, PIXELTYPE ** > traverser
Definition: basicimage.hxx:526
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read squared magnitude at offset from iterator position.
Definition: fftw3.hxx:1365
Definition: fftw3.hxx:1403
Accessor for items that are STL compatible vectors.
Definition: accessor.hxx:770
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
FFTWComplex & operator=(float o)
Definition: fftw3.hxx:291
Definition: basicimage.hxx:261
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read real part at offset from iterator position.
Definition: fftw3.hxx:1232
FFTWComplex< R > & operator+=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
add-assignment
Definition: fftw3.hxx:859
linalg::TemporaryMatrix< T > log10(MultiArrayView< 2, T, C > const &v)
FFTWComplex(FFTWComplex const &o)
Definition: fftw3.hxx:177
reference operator[](int i)
Definition: fftw3.hxx:373
value_type operator()(ITERATOR const &i) const
Read imaginary part at iterator position.
Definition: fftw3.hxx:1276
value_type const & const_reference
Definition: fftw3.hxx:148
bool operator!=(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
not equal
Definition: fftw3.hxx:841
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read imaginary part at offset from iterator position.
Definition: fftw3.hxx:1282
FFTWComplex(FFTWComplex< U > const &o)
Definition: fftw3.hxx:186
value_type operator()(ITERATOR const &i) const
Read phase at iterator position.
Definition: fftw3.hxx:1411
FFTWComplex & operator=(long double o)
Definition: fftw3.hxx:300
bool operator==(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
equal
Definition: fftw3.hxx:825
FFTWComplex & operator=(FFTWComplex< U > const &o)
Definition: fftw3.hxx:246
void combineTwoImages(...)
Combine two source images into destination image.
FFTWComplex(fftw_complex const &o)
Definition: fftw3.hxx:194
void copyImage(...)
Copy source image into destination image.
FFTWComplex< R > & operator*=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
multiply-assignment
Definition: fftw3.hxx:875
FFTWComplex(fftwl_complex const &o)
Definition: fftw3.hxx:210
value_type operator()(ITERATOR const &i) const
Read real part at iterator position.
Definition: fftw3.hxx:1226
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:939
IteratorTraits< const_traverser >::DefaultAccessor ConstAccessor
Definition: basicimage.hxx:576
value_type operator()(ITERATOR const &i) const
Read magnitude at iterator position.
Definition: fftw3.hxx:1385
value_type NormType
Definition: fftw3.hxx:160
Definition: fftw3.hxx:1268
ConstBasicImageIterator< PIXELTYPE, PIXELTYPE ** > const_traverser
Definition: basicimage.hxx:536
Definition: fftw3.hxx:1377
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
Fundamental class template for images.
Definition: basicimage.hxx:473
void set(value_type const &v, ITERATOR const &i) const
Write imaginary part at iterator position from a scalar.
Definition: fftw3.hxx:1288
FFTWComplex operator-() const
Definition: fftw3.hxx:353
void fourierTransformReal(...)
Real Fourier transforms for even and odd boundary conditions (aka. cosine and sine transforms)...
Definition: basicimage.hxx:293
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
linalg::TemporaryMatrix< T > tan(MultiArrayView< 2, T, C > const &v)
value_type * iterator
Definition: fftw3.hxx:152
const_reference operator[](int i) const
Definition: fftw3.hxx:378
FFTWComplex< R > & operator/=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
divide-assignment
Definition: fftw3.hxx:884
NormType magnitude() const
Definition: fftw3.hxx:363
LinearIntensityTransform< DestValueType, Multiplier > linearIntensityTransform(Multiplier scale, DestValueType offset)
Apply a linear transform to the source pixel values.
Definition: transformimage.hxx:681
void set(value_type const &v, ITERATOR const &i) const
Write real part at iterator position from a scalar.
Definition: fftw3.hxx:1238
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
T sign(T t)
Definition: mathutil.hxx:551
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1407
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read magnitude at offset from iterator position.
Definition: fftw3.hxx:1391
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
Wrapper class for the FFTW complex types 'fftw_complex'.
Definition: fftw3.hxx:131
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
value_type & reference
Definition: fftw3.hxx:144
FFTWComplex(std::complex< T > const &o)
Definition: fftw3.hxx:219
BasicImage< FFTWComplex<> > FFTWComplexImage
Definition: fftw3.hxx:1192
FFTWComplex(TinyVector< T, 2 > const &o)
Definition: fftw3.hxx:228
FFTWComplex & operator=(double o)
Definition: fftw3.hxx:282
FFTWComplex & operator=(fftwf_complex const &o)
Definition: fftw3.hxx:264
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.9.0 (Sun Aug 10 2014)