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

fftw3.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_FFTW3_HXX
37 #define VIGRA_FFTW3_HXX
38 
39 #include <cmath>
40 #include <functional>
41 #include <complex>
42 #include "stdimage.hxx"
43 #include "copyimage.hxx"
44 #include "transformimage.hxx"
45 #include "combineimages.hxx"
46 #include "numerictraits.hxx"
47 #include "imagecontainer.hxx"
48 #include <fftw3.h>
49 
50 namespace vigra {
51 
52 typedef double fftw_real;
53 
54 template <class T>
55 struct FFTWReal;
56 
57 template <>
58 struct FFTWReal<fftw_complex>
59 {
60  typedef double type;
61 };
62 
63 template <>
64 struct FFTWReal<fftwf_complex>
65 {
66  typedef float type;
67 };
68 
69 template <>
70 struct FFTWReal<fftwl_complex>
71 {
72  typedef long double type;
73 };
74 
75 template <class T>
76 struct FFTWReal2Complex;
77 
78 template <>
79 struct FFTWReal2Complex<double>
80 {
81  typedef fftw_complex type;
82  typedef fftw_plan plan_type;
83 };
84 
85 template <>
86 struct FFTWReal2Complex<float>
87 {
88  typedef fftwf_complex type;
89  typedef fftwf_plan plan_type;
90 };
91 
92 template <>
93 struct FFTWReal2Complex<long double>
94 {
95  typedef fftwl_complex type;
96  typedef fftwl_plan plan_type;
97 };
98 
99 /********************************************************/
100 /* */
101 /* FFTWComplex */
102 /* */
103 /********************************************************/
104 
105 /** \brief Wrapper class for the FFTW complex types '<TT>fftw_complex</TT>'.
106 
107  This class encapsulates the low-level complex number types provided by the
108  <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a> library (i.e.
109  '<TT>fftw_complex</TT>', '<TT>fftwf_complex</TT>', '<TT>fftwl_complex</TT>').
110  In particular, it provides constructors, member functions and
111  \ref FFTWComplexOperators "arithmetic operators" that make FFTW complex numbers
112  compatible with <tt>std::complex</tt>. In addition, the class defines
113  transformations to polar coordinates and \ref FFTWComplexAccessors "accessors".
114 
115  FFTWComplex implements the concepts \ref AlgebraicField and
116  \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt>
117  and <tt>FFTWComplexImage</tt> are defined.
118 
119  <b>See also:</b>
120  <ul>
121  <li> \ref FFTWComplexTraits
122  <li> \ref FFTWComplexOperators
123  <li> \ref FFTWComplexAccessors
124  </ul>
125 
126  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
127  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
128  Namespace: vigra
129 */
130 template <class Real = double>
132 {
133  public:
134  /** The wrapped complex type
135  */
136  typedef typename FFTWReal2Complex<Real>::type complex_type;
137 
138  /** The complex' component type, as defined in '<TT>fftw3.h</TT>'
139  */
140  typedef Real value_type;
141 
142  /** reference type (result of operator[])
143  */
144  typedef value_type & reference;
145 
146  /** const reference type (result of operator[] const)
147  */
148  typedef value_type const & const_reference;
149 
150  /** iterator type (result of begin() )
151  */
152  typedef value_type * iterator;
153 
154  /** const iterator type (result of begin() const)
155  */
156  typedef value_type const * const_iterator;
157 
158  /** The norm type (result of magnitude())
159  */
160  typedef value_type NormType;
161 
162  /** The squared norm type (result of squaredMagnitde())
163  */
164  typedef value_type SquaredNormType;
165 
166  /** Construct from real and imaginary part.
167  Default: 0.
168  */
169  FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0)
170  {
171  data_[0] = re;
172  data_[1] = im;
173  }
174 
175  /** Copy constructor.
176  */
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  */
358  SquaredNormType squaredMagnitude() const
359  { return data_[0]*data_[0]+data_[1]*data_[1]; }
360 
361  /** Magnitude (length of radius vector).
362  */
363  NormType magnitude() const
364  { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
365 
366  /** Phase angle.
367  */
368  value_type phase() const
369  { return VIGRA_CSTD::atan2(data_[1], data_[0]); }
370 
371  /** Access components as if number were a vector.
372  */
373  reference operator[](int i)
374  { return data_[i]; }
375 
376  /** Read components as if number were a vector.
377  */
378  const_reference operator[](int i) const
379  { return data_[i]; }
380 
381  /** Length of complex number (always 2).
382  */
383  int size() const
384  { return 2; }
385 
386  iterator begin()
387  { return data_; }
388 
389  iterator end()
390  { return data_ + 2; }
391 
392  const_iterator begin() const
393  { return data_; }
394 
395  const_iterator end() const
396  { return data_ + 2; }
397 
398  private:
399  complex_type data_;
400 };
401 
402 /********************************************************/
403 /* */
404 /* FFTWComplexTraits */
405 /* */
406 /********************************************************/
407 
408 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex
409 
410  The numeric and promote traits for fftw_complex and FFTWComplex follow
411  the general specifications for \ref NumericPromotionTraits and
412  \ref AlgebraicField. They are explicitly specialized for the types
413  involved:
414 
415  \code
416 
417  template<>
418  struct NumericTraits<fftw_complex>
419  {
420  typedef fftw_complex Promote;
421  typedef fftw_complex RealPromote;
422  typedef fftw_complex ComplexPromote;
423  typedef fftw_real ValueType;
424 
425  typedef VigraFalseType isIntegral;
426  typedef VigraFalseType isScalar;
427  typedef VigraFalseType isOrdered;
428  typedef VigraTrueType isComplex;
429 
430  // etc.
431  };
432 
433  template<class Real>
434  struct NumericTraits<FFTWComplex<Real> >
435  {
436  typedef FFTWComplex<Real> Promote;
437  typedef FFTWComplex<Real> RealPromote;
438  typedef FFTWComplex<Real> ComplexPromote;
439  typedef Real ValueType;
440 
441  typedef VigraFalseType isIntegral;
442  typedef VigraFalseType isScalar;
443  typedef VigraFalseType isOrdered;
444  typedef VigraTrueType isComplex;
445 
446  // etc.
447  };
448 
449  template<>
450  struct NormTraits<fftw_complex>
451  {
452  typedef fftw_complex Type;
453  typedef fftw_real SquaredNormType;
454  typedef fftw_real NormType;
455  };
456 
457  template<class Real>
458  struct NormTraits<FFTWComplex>
459  {
460  typedef FFTWComplex Type;
461  typedef fftw_real SquaredNormType;
462  typedef fftw_real NormType;
463  };
464 
465  template <>
466  struct PromoteTraits<fftw_complex, fftw_complex>
467  {
468  typedef fftw_complex Promote;
469  };
470 
471  template <>
472  struct PromoteTraits<fftw_complex, double>
473  {
474  typedef fftw_complex Promote;
475  };
476 
477  template <>
478  struct PromoteTraits<double, fftw_complex>
479  {
480  typedef fftw_complex Promote;
481  };
482 
483  template <class Real>
484  struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
485  {
486  typedef FFTWComplex<Real> Promote;
487  };
488 
489  template <class Real>
490  struct PromoteTraits<FFTWComplex<Real>, double>
491  {
492  typedef FFTWComplex<Real> Promote;
493  };
494 
495  template <class Real>
496  struct PromoteTraits<double, FFTWComplex<Real> >
497  {
498  typedef FFTWComplex<Real> Promote;
499  };
500  \endcode
501 
502  <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
503  <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
504  Namespace: vigra
505 
506 */
507 template<>
508 struct NumericTraits<fftw_complex>
509 {
510  typedef fftw_complex Type;
511  typedef fftw_complex Promote;
512  typedef fftw_complex RealPromote;
513  typedef fftw_complex ComplexPromote;
514  typedef fftw_real ValueType;
515 
516  typedef VigraFalseType isIntegral;
517  typedef VigraFalseType isScalar;
518  typedef NumericTraits<fftw_real>::isSigned isSigned;
519  typedef VigraFalseType isOrdered;
520  typedef VigraTrueType isComplex;
521 
522  static FFTWComplex<> zero() { return FFTWComplex<>(0.0, 0.0); }
523  static FFTWComplex<> one() { return FFTWComplex<>(1.0, 0.0); }
524  static FFTWComplex<> nonZero() { return one(); }
525 
526  static const Promote & toPromote(const Type & v) { return v; }
527  static const RealPromote & toRealPromote(const Type & v) { return v; }
528  static const Type & fromPromote(const Promote & v) { return v; }
529  static const Type & fromRealPromote(const RealPromote & v) { return v; }
530 };
531 
532 template<class Real>
533 struct NumericTraits<FFTWComplex<Real> >
534 {
535  typedef FFTWComplex<Real> Type;
536  typedef FFTWComplex<Real> Promote;
537  typedef FFTWComplex<Real> RealPromote;
538  typedef FFTWComplex<Real> ComplexPromote;
539  typedef typename Type::value_type ValueType;
540 
541  typedef VigraFalseType isIntegral;
542  typedef VigraFalseType isScalar;
543  typedef typename NumericTraits<ValueType>::isSigned isSigned;
544  typedef VigraFalseType isOrdered;
545  typedef VigraTrueType isComplex;
546 
547  static Type zero() { return Type(0.0, 0.0); }
548  static Type one() { return Type(1.0, 0.0); }
549  static Type nonZero() { return one(); }
550 
551  static const Promote & toPromote(const Type & v) { return v; }
552  static const RealPromote & toRealPromote(const Type & v) { return v; }
553  static const Type & fromPromote(const Promote & v) { return v; }
554  static const Type & fromRealPromote(const RealPromote & v) { return v; }
555 };
556 
557 template<>
558 struct NormTraits<fftw_complex>
559 {
560  typedef fftw_complex Type;
561  typedef fftw_real SquaredNormType;
562  typedef fftw_real NormType;
563 };
564 
565 template<class Real>
566 struct NormTraits<FFTWComplex<Real> >
567 {
568  typedef FFTWComplex<Real> Type;
569  typedef typename Type::SquaredNormType SquaredNormType;
570  typedef typename Type::NormType NormType;
571 };
572 
573 template <>
574 struct PromoteTraits<fftw_complex, fftw_complex>
575 {
576  typedef fftw_complex Promote;
577 };
578 
579 template <>
580 struct PromoteTraits<fftw_complex, double>
581 {
582  typedef fftw_complex Promote;
583 };
584 
585 template <>
586 struct PromoteTraits<double, fftw_complex>
587 {
588  typedef fftw_complex Promote;
589 };
590 
591 template <class Real>
592 struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
593 {
594  typedef FFTWComplex<Real> Promote;
595 };
596 
597 template <class Real>
598 struct PromoteTraits<FFTWComplex<Real>, double>
599 {
600  typedef FFTWComplex<Real> Promote;
601 };
602 
603 template <class Real>
604 struct PromoteTraits<double, FFTWComplex<Real> >
605 {
606  typedef FFTWComplex<Real> Promote;
607 };
608 
609 template<class T>
610 struct CanSkipInitialization<std::complex<T> >
611 {
612  typedef typename CanSkipInitialization<T>::type type;
613  static const bool value = type::asBool;
614 };
615 
616 template<class Real>
617 struct CanSkipInitialization<FFTWComplex<Real> >
618 {
619  typedef typename CanSkipInitialization<Real>::type type;
620  static const bool value = type::asBool;
621 };
622 
623 namespace multi_math {
624 
625 template <class ARG>
626 struct MultiMathOperand;
627 
628 template <class Real>
629 struct MultiMathOperand<FFTWComplex<Real> >
630 {
631  typedef MultiMathOperand<FFTWComplex<Real> > AllowOverload;
632  typedef FFTWComplex<Real> result_type;
633 
634  static const int ndim = 0;
635 
636  MultiMathOperand(FFTWComplex<Real> const & v)
637  : v_(v)
638  {}
639 
640  template <class SHAPE>
641  bool checkShape(SHAPE const &) const
642  {
643  return true;
644  }
645 
646  template <class SHAPE>
647  FFTWComplex<Real> const & operator[](SHAPE const &) const
648  {
649  return v_;
650  }
651 
652  void inc(unsigned int /*LEVEL*/) const
653  {}
654 
655  void reset(unsigned int /*LEVEL*/) const
656  {}
657 
658  FFTWComplex<Real> const & operator*() const
659  {
660  return v_;
661  }
662 
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>
966  a *= b;
967  return a;
968 }
969 
970  /// right multiplication with scalar double
971 template <class R>
973  a *= b;
974  return a;
975 }
976 
977  /// left multiplication with scalar double
978 template <class R>
980  b *= a;
981  return b;
982 }
983 
984  /// division
985 template <class R>
987  a /= b;
988  return a;
989 }
990 
991  /// right division with scalar double
992 template <class R>
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
1919 fourierTransform(FFTWComplexImage::const_traverser sul,
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
1960 fourierTransformInverse(FFTWComplexImage::const_traverser sul,
1963 {
1964  detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
1965 }
1966 
1967 template <class DestImageIterator, class DestAccessor>
1968 void fourierTransformInverse(FFTWComplexImage::const_traverser sul,
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
1983 fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
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 2D array views:
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  pass \ref ImageIterators and \ref DataAccessors :
2032  \code
2033  namespace vigra {
2034  template <class SrcImageIterator, class SrcAccessor,
2035  class FilterImageIterator, class FilterAccessor,
2036  class DestImageIterator, class DestAccessor>
2037  void applyFourierFilter(SrcImageIterator srcUpperLeft,
2038  SrcImageIterator srcLowerRight, SrcAccessor sa,
2039  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2040  DestImageIterator destUpperLeft, DestAccessor da);
2041  }
2042  \endcode
2043 
2044  use argument objects in conjunction with \ref ArgumentObjectFactories :
2045  \code
2046  namespace vigra {
2047  template <class SrcImageIterator, class SrcAccessor,
2048  class FilterImageIterator, class FilterAccessor,
2049  class DestImageIterator, class DestAccessor>
2050  void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2051  pair<FilterImageIterator, FilterAccessor> filter,
2052  pair<DestImageIterator, DestAccessor> dest);
2053  }
2054  \endcode
2055 
2056  <b> Usage:</b>
2057 
2058  <b>\#include</b> <vigra/fftw3.hxx><br>
2059  Namespace: vigra
2060 
2061  \code
2062  // create a Gaussian filter in Fourier space
2063  vigra::FImage gaussFilter(w, h), filter(w, h);
2064  for(int y=0; y<h; ++y)
2065  for(int x=0; x<w; ++x)
2066  {
2067  xx = float(x - w / 2) / w;
2068  yy = float(y - h / 2) / h;
2069 
2070  gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
2071  }
2072 
2073  // applyFourierFilter() expects the filter's DC in the upper left
2074  moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
2075 
2076  vigra::FFTWComplexImage result(w, h);
2077 
2078  vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
2079  \endcode
2080 
2081  For inspection of the result, \ref FFTWMagnitudeAccessor might be
2082  useful. If you want to apply the same filter repeatedly, it may be more
2083  efficient to use the FFTW functions directly with FFTW plans optimized
2084  for good performance.
2085 */
2086 doxygen_overloaded_function(template <...> void applyFourierFilter)
2087 
2088 template <class SrcImageIterator, class SrcAccessor,
2089  class FilterImageIterator, class FilterAccessor,
2090  class DestImageIterator, class DestAccessor>
2091 void applyFourierFilter(SrcImageIterator srcUpperLeft,
2092  SrcImageIterator srcLowerRight, SrcAccessor sa,
2093  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2094  DestImageIterator destUpperLeft, DestAccessor da)
2095 {
2096  // copy real input images into a complex one...
2097  int w = int(srcLowerRight.x - srcUpperLeft.x);
2098  int h = int(srcLowerRight.y - srcUpperLeft.y);
2099 
2100  FFTWComplexImage workImage(w, h);
2101  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2102  destImage(workImage, FFTWWriteRealAccessor<>()));
2103 
2104  // ...and call the impl
2105  FFTWComplexImage const & cworkImage = workImage;
2106  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2107  filterUpperLeft, fa,
2108  destUpperLeft, da);
2109 }
2110 
2111 template <class FilterImageIterator, class FilterAccessor,
2112  class DestImageIterator, class DestAccessor>
2113 inline
2114 void applyFourierFilter(
2115  FFTWComplexImage::const_traverser srcUpperLeft,
2116  FFTWComplexImage::const_traverser srcLowerRight,
2118  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2119  DestImageIterator destUpperLeft, DestAccessor da)
2120 {
2121  int w = srcLowerRight.x - srcUpperLeft.x;
2122  int h = srcLowerRight.y - srcUpperLeft.y;
2123 
2124  // test for right memory layout (fftw expects a 2*width*height floats array)
2125  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2126  applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
2127  filterUpperLeft, fa,
2128  destUpperLeft, da);
2129  else
2130  {
2131  FFTWComplexImage workImage(w, h);
2132  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2133  destImage(workImage));
2134 
2135  FFTWComplexImage const & cworkImage = workImage;
2136  applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2137  filterUpperLeft, fa,
2138  destUpperLeft, da);
2139  }
2140 }
2141 
2142 template <class SrcImageIterator, class SrcAccessor,
2143  class FilterImageIterator, class FilterAccessor,
2144  class DestImageIterator, class DestAccessor>
2145 inline
2146 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2147  pair<FilterImageIterator, FilterAccessor> filter,
2148  pair<DestImageIterator, DestAccessor> dest)
2149 {
2150  applyFourierFilter(src.first, src.second, src.third,
2151  filter.first, filter.second,
2152  dest.first, dest.second);
2153 }
2154 
2155 template <class FilterImageIterator, class FilterAccessor,
2156  class DestImageIterator, class DestAccessor>
2157 void applyFourierFilterImpl(
2158  FFTWComplexImage::const_traverser srcUpperLeft,
2159  FFTWComplexImage::const_traverser srcLowerRight,
2161  FilterImageIterator filterUpperLeft, FilterAccessor fa,
2162  DestImageIterator destUpperLeft, DestAccessor da)
2163 {
2164  int w = int(srcLowerRight.x - srcUpperLeft.x);
2165  int h = int(srcLowerRight.y - srcUpperLeft.y);
2166 
2167  FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
2168 
2169  // FFT from srcImage to complexResultImg
2170  fftw_plan forwardPlan=
2171  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2172  (fftw_complex *)complexResultImg.begin(),
2173  FFTW_FORWARD, FFTW_ESTIMATE );
2174  fftw_execute(forwardPlan);
2175  fftw_destroy_plan(forwardPlan);
2176 
2177  // convolve in freq. domain (in complexResultImg)
2178  combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
2179  destImage(complexResultImg), std::multiplies<FFTWComplex<> >());
2180 
2181  // FFT back into spatial domain (inplace in complexResultImg)
2182  fftw_plan backwardPlan=
2183  fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
2184  (fftw_complex *)complexResultImg.begin(),
2185  FFTW_BACKWARD, FFTW_ESTIMATE);
2186  fftw_execute(backwardPlan);
2187  fftw_destroy_plan(backwardPlan);
2188 
2189  typedef typename
2190  NumericTraits<typename DestAccessor::value_type>::isScalar
2191  isScalarResult;
2192 
2193  // normalization (after FFTs), maybe stripping imaginary part
2194  applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
2195  isScalarResult());
2196 }
2197 
2198 template <class DestImageIterator, class DestAccessor>
2199 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
2200  DestImageIterator destUpperLeft,
2201  DestAccessor da,
2202  VigraFalseType)
2203 {
2204  double normFactor= 1.0/(srcImage.width() * srcImage.height());
2205 
2206  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2207  {
2208  DestImageIterator dIt= destUpperLeft;
2209  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
2210  {
2211  da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
2212  da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
2213  }
2214  }
2215 }
2216 
2217 inline
2218 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
2219  FFTWComplexImage::traverser destUpperLeft,
2221  VigraFalseType)
2222 {
2223  transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
2224  linearIntensityTransform<FFTWComplex<> >(1.0/(srcImage.width() * srcImage.height())));
2225 }
2226 
2227 template <class DestImageIterator, class DestAccessor>
2228 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
2229  DestImageIterator destUpperLeft,
2230  DestAccessor da,
2231  VigraTrueType)
2232 {
2233  double normFactor= 1.0/(srcImage.width() * srcImage.height());
2234 
2235  for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2236  {
2237  DestImageIterator dIt= destUpperLeft;
2238  for(int x= 0; x< srcImage.width(); x++, dIt.x++)
2239  da.set(srcImage(x, y).re()*normFactor, dIt);
2240  }
2241 }
2242 
2243 /**********************************************************/
2244 /* */
2245 /* applyFourierFilterFamily */
2246 /* */
2247 /**********************************************************/
2248 
2249 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
2250 
2251  This provides the same functionality as \ref applyFourierFilter(),
2252  but applying several filters at once allows to avoid
2253  repeated Fourier transforms of the source image.
2254 
2255  Filters and result images must be stored in \ref vigra::ImageArray data
2256  structures. In contrast to \ref applyFourierFilter(), this function adjusts
2257  the size of the result images and the the length of the array.
2258 
2259  <b> Declarations:</b>
2260 
2261  pass 2D array views:
2262  \code
2263  namespace vigra {
2264  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2265  void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2266  SrcImageIterator srcLowerRight, SrcAccessor sa,
2267  const ImageArray<FilterType> &filters,
2268  ImageArray<FFTWComplexImage> &results)
2269  }
2270  \endcode
2271 
2272  pass \ref ImageIterators and \ref DataAccessors :
2273  \code
2274  namespace vigra {
2275  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2276  void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2277  SrcImageIterator srcLowerRight, SrcAccessor sa,
2278  const ImageArray<FilterType> &filters,
2279  ImageArray<FFTWComplexImage> &results)
2280  }
2281  \endcode
2282 
2283  use argument objects in conjunction with \ref ArgumentObjectFactories :
2284  \code
2285  namespace vigra {
2286  template <class SrcImageIterator, class SrcAccessor, class FilterType>
2287  void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2288  const ImageArray<FilterType> &filters,
2289  ImageArray<FFTWComplexImage> &results)
2290  }
2291  \endcode
2292 
2293  <b> Usage:</b>
2294 
2295  <b>\#include</b> <vigra/fftw3.hxx><br>
2296  Namespace: vigra
2297 
2298  \code
2299  // assuming the presence of a real-valued image named "image" to
2300  // be filtered in this example
2301 
2302  vigra::ImageArray<vigra::FImage> filters(16, image.size());
2303 
2304  for (int i=0; i<filters.size(); i++)
2305  // create some meaningful filters here
2306  createMyFilterOfScale(i, destImage(filters[i]));
2307 
2308  vigra::ImageArray<vigra::FFTWComplexImage> results();
2309 
2310  vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
2311  \endcode
2312 */
2313 doxygen_overloaded_function(template <...> void applyFourierFilterFamily)
2314 
2315 template <class SrcImageIterator, class SrcAccessor,
2316  class FilterType, class DestImage>
2317 inline
2318 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2319  const ImageArray<FilterType> &filters,
2320  ImageArray<DestImage> &results)
2321 {
2322  applyFourierFilterFamily(src.first, src.second, src.third,
2323  filters, results);
2324 }
2325 
2326 template <class SrcImageIterator, class SrcAccessor,
2327  class FilterType, class DestImage>
2328 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2329  SrcImageIterator srcLowerRight, SrcAccessor sa,
2330  const ImageArray<FilterType> &filters,
2331  ImageArray<DestImage> &results)
2332 {
2333  int w = int(srcLowerRight.x - srcUpperLeft.x);
2334  int h = int(srcLowerRight.y - srcUpperLeft.y);
2335 
2336  FFTWComplexImage workImage(w, h);
2337  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2338  destImage(workImage, FFTWWriteRealAccessor<>()));
2339 
2340  FFTWComplexImage const & cworkImage = workImage;
2341  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2342  filters, results);
2343 }
2344 
2345 template <class FilterType, class DestImage>
2346 inline
2347 void applyFourierFilterFamily(
2348  FFTWComplexImage::const_traverser srcUpperLeft,
2349  FFTWComplexImage::const_traverser srcLowerRight,
2351  const ImageArray<FilterType> &filters,
2352  ImageArray<DestImage> &results)
2353 {
2354  int w= srcLowerRight.x - srcUpperLeft.x;
2355 
2356  // test for right memory layout (fftw expects a 2*width*height floats array)
2357  if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2358  applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
2359  filters, results);
2360  else
2361  {
2362  int h = srcLowerRight.y - srcUpperLeft.y;
2363  FFTWComplexImage workImage(w, h);
2364  copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2365  destImage(workImage));
2366 
2367  FFTWComplexImage const & cworkImage = workImage;
2368  applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2369  filters, results);
2370  }
2371 }
2372 
2373 template <class FilterType, class DestImage>
2374 void applyFourierFilterFamilyImpl(
2375  FFTWComplexImage::const_traverser srcUpperLeft,
2376  FFTWComplexImage::const_traverser srcLowerRight,
2378  const ImageArray<FilterType> &filters,
2379  ImageArray<DestImage> &results)
2380 {
2381  // FIXME: sa is not used
2382  // (maybe check if StandardAccessor, else copy?)
2383 
2384  // make sure the filter images have the right dimensions
2385  vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
2386  "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2387 
2388  // make sure the result image array has the right dimensions
2389  results.resize(filters.size());
2390  results.resizeImages(filters.imageSize());
2391 
2392  // FFT from srcImage to freqImage
2393  int w = int(srcLowerRight.x - srcUpperLeft.x);
2394  int h = int(srcLowerRight.y - srcUpperLeft.y);
2395 
2396  FFTWComplexImage freqImage(w, h);
2397  FFTWComplexImage result(w, h);
2398 
2399  fftw_plan forwardPlan=
2400  fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2401  (fftw_complex *)freqImage.begin(),
2402  FFTW_FORWARD, FFTW_ESTIMATE );
2403  fftw_execute(forwardPlan);
2404  fftw_destroy_plan(forwardPlan);
2405 
2406  fftw_plan backwardPlan=
2407  fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
2408  (fftw_complex *)result.begin(),
2409  FFTW_BACKWARD, FFTW_ESTIMATE );
2410  typedef typename
2411  NumericTraits<typename DestImage::Accessor::value_type>::isScalar
2412  isScalarResult;
2413 
2414  // convolve with filters in freq. domain
2415  for (unsigned int i= 0; i < filters.size(); i++)
2416  {
2417  combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
2418  destImage(result), std::multiplies<FFTWComplex<> >());
2419 
2420  // FFT back into spatial domain (inplace in destImage)
2421  fftw_execute(backwardPlan);
2422 
2423  // normalization (after FFTs), maybe stripping imaginary part
2424  applyFourierFilterImplNormalization(result,
2425  results[i].upperLeft(), results[i].accessor(),
2426  isScalarResult());
2427  }
2428  fftw_destroy_plan(backwardPlan);
2429 }
2430 
2431 /********************************************************/
2432 /* */
2433 /* fourierTransformReal */
2434 /* */
2435 /********************************************************/
2436 
2437 /** \brief Real Fourier transforms for even and odd boundary conditions
2438  (aka. cosine and sine transforms).
2439 
2440 
2441  If the image is real and has even symmetry, its Fourier transform
2442  is also real and has even symmetry. The Fourier transform of a real image with odd
2443  symmetry is imaginary and has odd symmetry. In either case, only about a quarter
2444  of the pixels need to be stored because the rest can be calculated from the symmetry
2445  properties. This is especially useful, if the original image is implicitly assumed
2446  to have reflective or anti-reflective boundary conditions. Then the "negative"
2447  pixel locations are defined as
2448 
2449  \code
2450  even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1)
2451  odd (anti-reflective boundary conditions): f[-1] = 0
2452  f[-x] = -f[x-2] (x = 2,...,N-1)
2453  \endcode
2454 
2455  end similar at the other boundary (see the FFTW documentation for details).
2456  This has the advantage that more efficient Fourier transforms that use only
2457  real numbers can be implemented. These are also known as cosine and sine transforms
2458  respectively.
2459 
2460  If you use the odd transform it is important to note that in the Fourier domain,
2461  the DC component is always zero and is therefore dropped from the data structure.
2462  This means that index 0 in an odd symmetric Fourier domain image refers to
2463  the <i>first</i> harmonic. This is especially important if an image is first
2464  cosine transformed (even symmetry), then in the Fourier domain multiplied
2465  with an odd symmetric filter (e.g. a first derivative) and finally transformed
2466  back to the spatial domain with a sine transform (odd symmetric). For this to work
2467  properly the image must be shifted left or up by one pixel (depending on whether
2468  the x- or y-axis is odd symmetric) before the inverse transform can be applied.
2469  (see example below).
2470 
2471  The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
2472  where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
2473  whether the x- and y-axis is to be transformed using even or odd symmetry.
2474  The same functions can be used for both the forward and inverse transforms,
2475  only the normalization changes. For signal processing, the following
2476  normalization factors are most appropriate:
2477 
2478  \code
2479  forward inverse
2480  ------------------------------------------------------------
2481  X even, Y even 1.0 4.0 * (w-1) * (h-1)
2482  X even, Y odd -1.0 -4.0 * (w-1) * (h+1)
2483  X odd, Y even -1.0 -4.0 * (w+1) * (h-1)
2484  X odd, Y odd 1.0 4.0 * (w+1) * (h+1)
2485  \endcode
2486 
2487  where <tt>w</tt> and <tt>h</tt> denote the image width and height.
2488 
2489  <b> Declarations:</b>
2490 
2491  pass 2D array views:
2492  \code
2493  namespace vigra {
2494  template <class SrcTraverser, class SrcAccessor,
2495  class DestTraverser, class DestAccessor>
2496  void
2497  fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2498  DestTraverser dul, DestAccessor dest, fftw_real norm);
2499 
2500  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2501  }
2502  \endcode
2503 
2504  pass \ref ImageIterators and \ref DataAccessors :
2505  \code
2506  namespace vigra {
2507  template <class SrcTraverser, class SrcAccessor,
2508  class DestTraverser, class DestAccessor>
2509  void
2510  fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2511  DestTraverser dul, DestAccessor dest, fftw_real norm);
2512 
2513  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2514  }
2515  \endcode
2516 
2517 
2518  use argument objects in conjunction with \ref ArgumentObjectFactories :
2519  \code
2520  namespace vigra {
2521  template <class SrcTraverser, class SrcAccessor,
2522  class DestTraverser, class DestAccessor>
2523  void
2524  fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2525  pair<DestTraverser, DestAccessor> dest, fftw_real norm);
2526 
2527  fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
2528  }
2529  \endcode
2530 
2531  <b> Usage:</b>
2532 
2533  <b>\#include</b> <vigra/fftw3.hxx><br>
2534  Namespace: vigra
2535 
2536  \code
2537  vigra::FImage spatial(width,height), fourier(width,height);
2538  ... // fill image with data
2539 
2540  // forward cosine transform == reflective boundary conditions
2541  fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
2542 
2543  // multiply with a first derivative of Gaussian in x-direction
2544  for(int y = 0; y < height; ++y)
2545  {
2546  for(int x = 1; x < width; ++x)
2547  {
2548  double dx = x * M_PI / (width - 1);
2549  double dy = y * M_PI / (height - 1);
2550  fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
2551  }
2552  fourier(width-1, y) = 0.0;
2553  }
2554 
2555  // inverse transform -- odd symmetry in x-direction, even in y,
2556  // due to symmetry of the filter
2557  fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
2558  (fftw_real)-4.0 * (width+1) * (height-1));
2559  \endcode
2560 */
2561 doxygen_overloaded_function(template <...> void fourierTransformReal)
2562 
2563 template <class SrcTraverser, class SrcAccessor,
2564  class DestTraverser, class DestAccessor>
2565 inline void
2566 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2567  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2568 {
2569  fourierTransformRealEE(src.first, src.second, src.third,
2570  dest.first, dest.second, norm);
2571 }
2572 
2573 template <class SrcTraverser, class SrcAccessor,
2574  class DestTraverser, class DestAccessor>
2575 inline void
2576 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2577  DestTraverser dul, DestAccessor dest, fftw_real norm)
2578 {
2579  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2580  norm, FFTW_REDFT00, FFTW_REDFT00);
2581 }
2582 
2583 template <class DestTraverser, class DestAccessor>
2584 inline void
2585 fourierTransformRealEE(
2589  DestTraverser dul, DestAccessor dest, fftw_real norm)
2590 {
2591  int w = slr.x - sul.x;
2592 
2593  // test for right memory layout (fftw expects a width*height fftw_real array)
2594  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2595  fourierTransformRealImpl(sul, slr, dul, dest,
2596  norm, FFTW_REDFT00, FFTW_REDFT00);
2597  else
2598  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2599  norm, FFTW_REDFT00, FFTW_REDFT00);
2600 }
2601 
2602 /********************************************************************/
2603 
2604 template <class SrcTraverser, class SrcAccessor,
2605  class DestTraverser, class DestAccessor>
2606 inline void
2607 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2608  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2609 {
2610  fourierTransformRealOE(src.first, src.second, src.third,
2611  dest.first, dest.second, norm);
2612 }
2613 
2614 template <class SrcTraverser, class SrcAccessor,
2615  class DestTraverser, class DestAccessor>
2616 inline void
2617 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2618  DestTraverser dul, DestAccessor dest, fftw_real norm)
2619 {
2620  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2621  norm, FFTW_RODFT00, FFTW_REDFT00);
2622 }
2623 
2624 template <class DestTraverser, class DestAccessor>
2625 inline void
2626 fourierTransformRealOE(
2630  DestTraverser dul, DestAccessor dest, fftw_real norm)
2631 {
2632  int w = slr.x - sul.x;
2633 
2634  // test for right memory layout (fftw expects a width*height fftw_real array)
2635  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2636  fourierTransformRealImpl(sul, slr, dul, dest,
2637  norm, FFTW_RODFT00, FFTW_REDFT00);
2638  else
2639  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2640  norm, FFTW_RODFT00, FFTW_REDFT00);
2641 }
2642 
2643 /********************************************************************/
2644 
2645 template <class SrcTraverser, class SrcAccessor,
2646  class DestTraverser, class DestAccessor>
2647 inline void
2648 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2649  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2650 {
2651  fourierTransformRealEO(src.first, src.second, src.third,
2652  dest.first, dest.second, norm);
2653 }
2654 
2655 template <class SrcTraverser, class SrcAccessor,
2656  class DestTraverser, class DestAccessor>
2657 inline void
2658 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2659  DestTraverser dul, DestAccessor dest, fftw_real norm)
2660 {
2661  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2662  norm, FFTW_REDFT00, FFTW_RODFT00);
2663 }
2664 
2665 template <class DestTraverser, class DestAccessor>
2666 inline void
2667 fourierTransformRealEO(
2671  DestTraverser dul, DestAccessor dest, fftw_real norm)
2672 {
2673  int w = slr.x - sul.x;
2674 
2675  // test for right memory layout (fftw expects a width*height fftw_real array)
2676  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2677  fourierTransformRealImpl(sul, slr, dul, dest,
2678  norm, FFTW_REDFT00, FFTW_RODFT00);
2679  else
2680  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2681  norm, FFTW_REDFT00, FFTW_RODFT00);
2682 }
2683 
2684 /********************************************************************/
2685 
2686 template <class SrcTraverser, class SrcAccessor,
2687  class DestTraverser, class DestAccessor>
2688 inline void
2689 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2690  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2691 {
2692  fourierTransformRealOO(src.first, src.second, src.third,
2693  dest.first, dest.second, norm);
2694 }
2695 
2696 template <class SrcTraverser, class SrcAccessor,
2697  class DestTraverser, class DestAccessor>
2698 inline void
2699 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2700  DestTraverser dul, DestAccessor dest, fftw_real norm)
2701 {
2702  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2703  norm, FFTW_RODFT00, FFTW_RODFT00);
2704 }
2705 
2706 template <class DestTraverser, class DestAccessor>
2707 inline void
2708 fourierTransformRealOO(
2712  DestTraverser dul, DestAccessor dest, fftw_real norm)
2713 {
2714  int w = slr.x - sul.x;
2715 
2716  // test for right memory layout (fftw expects a width*height fftw_real array)
2717  if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2718  fourierTransformRealImpl(sul, slr, dul, dest,
2719  norm, FFTW_RODFT00, FFTW_RODFT00);
2720  else
2721  fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2722  norm, FFTW_RODFT00, FFTW_RODFT00);
2723 }
2724 
2725 /*******************************************************************/
2726 
2727 template <class SrcTraverser, class SrcAccessor,
2728  class DestTraverser, class DestAccessor>
2729 void
2730 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2731  DestTraverser dul, DestAccessor dest,
2732  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2733 {
2734  FFTWRealImage workImage(slr - sul);
2735  copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2736  FFTWRealImage const & cworkImage = workImage;
2737  fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
2738  dul, dest, norm, kindx, kindy);
2739 }
2740 
2741 
2742 template <class DestTraverser, class DestAccessor>
2743 void
2744 fourierTransformRealImpl(
2747  DestTraverser dul, DestAccessor dest,
2748  fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2749 {
2750  int w = slr.x - sul.x;
2751  int h = slr.y - sul.y;
2752  BasicImage<fftw_real> res(w, h);
2753 
2754  fftw_plan plan = fftw_plan_r2r_2d(h, w,
2755  (fftw_real *)&(*sul), (fftw_real *)res.begin(),
2756  kindy, kindx, FFTW_ESTIMATE);
2757  fftw_execute(plan);
2758  fftw_destroy_plan(plan);
2759 
2760  if(norm != 1.0)
2761  transformImage(srcImageRange(res), destIter(dul, dest),
2762  std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2763  else
2764  copyImage(srcImageRange(res), destIter(dul, dest));
2765 }
2766 
2767 
2768 //@}
2769 
2770 } // namespace vigra
2771 
2772 #endif // VIGRA_FFTW3_HXX
FFTWComplex(fftwf_complex const &o)
Definition: fftw3.hxx:202
Size2D imageSize() const
Definition: imagecontainer.hxx:411
iterator begin()
Definition: basicimage.hxx:963
Definition: fftw3.hxx:1217
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&#39;s upper left...
FFTWComplex< R > & operator*=(FFTWComplex< R > &a, double b)
multiply-assignment with scalar double
Definition: fftw3.hxx:908
FFTWComplex< R > & operator/=(FFTWComplex< R > &a, double b)
divide-assignment with scalar double
Definition: fftw3.hxx:916
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1654
BasicImage< fftw_real > FFTWRealImage
Definition: fftw3.hxx:1132
Fundamental class template for arrays of equal-sized images.
Definition: imagecontainer.hxx:72
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
value_type SquaredNormType
Definition: fftw3.hxx:164
Accessor accessor()
Definition: basicimage.hxx:1064
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
Real value_type
The accessor&#39;s value type.
Definition: fftw3.hxx:1381
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
Real value_type
The accessor&#39;s value type.
Definition: fftw3.hxx:1324
Definition: fftw3.hxx:1319
value_type phase() const
Definition: fftw3.hxx:368
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.
int width() const
Definition: basicimage.hxx:838
Real value_type
The accessor&#39;s value type.
Definition: fftw3.hxx:1355
Definition: array_vector.hxx:903
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
FFTWComplex< R > & operator+=(FFTWComplex< R > &a, double b)
add-assignment with scalar double
Definition: fftw3.hxx:894
Two dimensional difference vector.
Definition: diff2d.hxx:185
Real value_type
The accessor&#39;s value type.
Definition: fftw3.hxx:1222
FFTWComplex & operator=(std::complex< T > const &o)
Definition: fftw3.hxx:320
Real value_type
The accessor&#39;s value type.
Definition: fftw3.hxx:1272
FFTWComplex & operator=(FFTWComplex const &o)
Definition: fftw3.hxx:236
Definition: accessor.hxx:43
value_type operator()(ITERATOR const &i) const
Read squared magnitude at iterator position.
Definition: fftw3.hxx:1359
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
int height() const
Definition: basicimage.hxx:845
Accessor for items that are STL compatible vectors.
Definition: accessor.hxx:771
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
FFTWComplex< R > & operator-=(FFTWComplex< R > &a, double b)
subtract-assignment with scalar double
Definition: fftw3.hxx:901
FFTWComplex & operator=(float o)
Definition: fftw3.hxx:291
Definition: basicimage.hxx:262
bool operator!=(double a, const FFTWComplex< R > &b)
not equal
Definition: fftw3.hxx:853
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read real part at offset from iterator position.
Definition: fftw3.hxx:1232
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
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
void resize(int width, int height)
Definition: basicimage.hxx:776
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
size_type size() const
Definition: imagecontainer.hxx:228
void copyImage(...)
Copy source image into destination image.
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
FFTWComplex< R > operator*(double a, FFTWComplex< R > b)
left multiplication with scalar double
Definition: fftw3.hxx:979
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
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
traverser lowerRight()
Definition: basicimage.hxx:934
Definition: fftw3.hxx:1377
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
Fundamental class template for images.
Definition: basicimage.hxx:473
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:294
FFTWComplex< R > operator-(double a, FFTWComplex< R > const &b)
left subtraction with scalar double
Definition: fftw3.hxx:959
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
void resize(size_type newSize)
Definition: imagecontainer.hxx:310
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
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:800
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
T sign(T t)
The sign function.
Definition: mathutil.hxx:553
virtual void resizeImages(const Diff2D &newSize)
Definition: imagecontainer.hxx:418
Real value_type
The accessor&#39;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
Wrapper class for the FFTW complex types &#39;fftw_complex&#39;.
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< R > operator/(FFTWComplex< R > a, double b)
right division with scalar double
Definition: fftw3.hxx:993
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
traverser upperLeft()
Definition: basicimage.hxx:923
FFTWComplex< R > operator+(double a, FFTWComplex< R > b)
left addition with scalar double
Definition: fftw3.hxx:938
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.10.0