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

vigra/fftw3.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_FFTW3_HXX
00037 #define VIGRA_FFTW3_HXX
00038 
00039 #include <cmath>
00040 #include <functional>
00041 #include <complex>
00042 #include "stdimage.hxx"
00043 #include "copyimage.hxx"
00044 #include "transformimage.hxx"
00045 #include "combineimages.hxx"
00046 #include "numerictraits.hxx"
00047 #include "imagecontainer.hxx"
00048 #include <fftw3.h>
00049 
00050 namespace vigra {
00051 
00052 typedef double fftw_real;
00053 
00054 template <class T>
00055 struct FFTWReal;
00056 
00057 template <>
00058 struct FFTWReal<fftw_complex>
00059 {
00060     typedef double type;
00061 };
00062 
00063 template <>
00064 struct FFTWReal<fftwf_complex>
00065 {
00066     typedef float type;
00067 };
00068 
00069 template <>
00070 struct FFTWReal<fftwl_complex>
00071 {
00072     typedef long double type;
00073 };
00074 
00075 template <class T>
00076 struct FFTWReal2Complex;
00077 
00078 template <>
00079 struct FFTWReal2Complex<double>
00080 {
00081     typedef fftw_complex type;
00082     typedef fftw_plan plan_type;
00083 };
00084 
00085 template <>
00086 struct FFTWReal2Complex<float>
00087 {
00088     typedef fftwf_complex type;
00089     typedef fftwf_plan plan_type;
00090 };
00091 
00092 template <>
00093 struct FFTWReal2Complex<long double>
00094 {
00095     typedef fftwl_complex type;
00096     typedef fftwl_plan plan_type;
00097 };
00098 
00099 /********************************************************/
00100 /*                                                      */
00101 /*                    FFTWComplex                       */
00102 /*                                                      */
00103 /********************************************************/
00104 
00105 /** \brief Wrapper class for the FFTW complex types '<TT>fftw_complex</TT>'.
00106 
00107     This class encapsulates the low-level complex number types provided by the 
00108     <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a> library (i.e. 
00109     '<TT>fftw_complex</TT>', '<TT>fftwf_complex</TT>', '<TT>fftwl_complex</TT>'). 
00110     In particular, it provides constructors, member functions and 
00111     \ref FFTWComplexOperators "arithmetic operators" that make FFTW complex numbers
00112     compatible with <tt>std::complex</tt>. In addition, the class defines 
00113     transformations to polar coordinates and \ref FFTWComplexAccessors "accessors".
00114 
00115     FFTWComplex implements the concepts \ref AlgebraicField and
00116     \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt>
00117     and <tt>FFTWComplexImage</tt> are defined.
00118 
00119     <b>See also:</b>
00120     <ul>
00121         <li> \ref FFTWComplexTraits
00122         <li> \ref FFTWComplexOperators
00123         <li> \ref FFTWComplexAccessors
00124     </ul>
00125 
00126     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
00127     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
00128     Namespace: vigra
00129 */
00130 template <class Real = double>
00131 class FFTWComplex
00132 {
00133   public:
00134         /** The wrapped complex type
00135         */
00136       typedef typename FFTWReal2Complex<Real>::type complex_type;
00137 
00138         /** The complex' component type, as defined in '<TT>fftw3.h</TT>'
00139         */
00140     typedef Real value_type;
00141 
00142         /** reference type (result of operator[])
00143         */
00144     typedef value_type & reference;
00145 
00146         /** const reference type (result of operator[] const)
00147         */
00148     typedef value_type const & const_reference;
00149 
00150         /** iterator type (result of begin() )
00151         */
00152     typedef value_type * iterator;
00153 
00154         /** const iterator type (result of begin() const)
00155         */
00156     typedef value_type const * const_iterator;
00157 
00158         /** The norm type (result of magnitude())
00159         */
00160     typedef value_type NormType;
00161 
00162         /** The squared norm type (result of squaredMagnitde())
00163         */
00164     typedef value_type SquaredNormType;
00165 
00166         /** Construct from real and imaginary part.
00167             Default: 0.
00168         */
00169     FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0)
00170     {
00171         data_[0] = re;
00172         data_[1] = im;
00173     }
00174 
00175         /** Copy constructor.
00176         */
00177     FFTWComplex(FFTWComplex const & o)
00178     {
00179         data_[0] = o.data_[0];
00180         data_[1] = o.data_[1];
00181     }
00182 
00183         /** Copy constructor.
00184         */
00185     template <class U>
00186     FFTWComplex(FFTWComplex<U> const & o)
00187     {
00188         data_[0] = (Real)o.real();
00189         data_[1] = (Real)o.imag();
00190     }
00191 
00192         /** Construct from plain <TT>fftw_complex</TT>.
00193         */
00194     FFTWComplex(fftw_complex const & o)
00195     {
00196         data_[0] = (Real)o[0];
00197         data_[1] = (Real)o[1];
00198     }
00199 
00200         /** Construct from plain <TT>fftwf_complex</TT>.
00201         */
00202     FFTWComplex(fftwf_complex const & o)
00203     {
00204         data_[0] = (Real)o[0];
00205         data_[1] = (Real)o[1];
00206     }
00207 
00208         /** Construct from plain <TT>fftwl_complex</TT>.
00209         */
00210     FFTWComplex(fftwl_complex const & o)
00211     {
00212         data_[0] = (Real)o[0];
00213         data_[1] = (Real)o[1];
00214     }
00215 
00216         /** Construct from std::complex.
00217         */
00218     template <class T>
00219     FFTWComplex(std::complex<T> const & o)
00220     {
00221         data_[0] = (Real)o.real();
00222         data_[1] = (Real)o.imag();
00223     }
00224 
00225         /** Construct from TinyVector.
00226         */
00227     template <class T>
00228     FFTWComplex(TinyVector<T, 2> const & o)
00229     {
00230         data_[0] = (Real)o[0];
00231         data_[1] = (Real)o[1];
00232     }
00233 
00234         /** Assignment.
00235         */
00236     FFTWComplex& operator=(FFTWComplex const & o)
00237     {
00238         data_[0] = o.data_[0];
00239         data_[1] = o.data_[1];
00240         return *this;
00241     }
00242 
00243         /** Assignment.
00244         */
00245     template <class U>
00246     FFTWComplex& operator=(FFTWComplex<U> const & o)
00247     {
00248         data_[0] = (Real)o.real();
00249         data_[1] = (Real)o.imag();
00250         return *this;
00251     }
00252 
00253         /** Assignment.
00254         */
00255     FFTWComplex& operator=(fftw_complex const & o)
00256     {
00257         data_[0] = (Real)o[0];
00258         data_[1] = (Real)o[1];
00259         return *this;
00260     }
00261 
00262         /** Assignment.
00263         */
00264     FFTWComplex& operator=(fftwf_complex const & o)
00265     {
00266         data_[0] = (Real)o[0];
00267         data_[1] = (Real)o[1];
00268         return *this;
00269     }
00270 
00271         /** Assignment.
00272         */
00273     FFTWComplex& operator=(fftwl_complex const & o)
00274     {
00275         data_[0] = (Real)o[0];
00276         data_[1] = (Real)o[1];
00277         return *this;
00278     }
00279 
00280         /** Assignment.
00281         */
00282     FFTWComplex& operator=(double o)
00283     {
00284         data_[0] = (Real)o;
00285         data_[1] = 0.0;
00286         return *this;
00287     }
00288 
00289         /** Assignment.
00290         */
00291     FFTWComplex& operator=(float o)
00292     {
00293         data_[0] = (Real)o;
00294         data_[1] = 0.0;
00295         return *this;
00296     }
00297 
00298         /** Assignment.
00299         */
00300     FFTWComplex& operator=(long double o)
00301     {
00302         data_[0] = (Real)o;
00303         data_[1] = 0.0;
00304         return *this;
00305     }
00306 
00307         /** Assignment.
00308         */
00309     template <class T>
00310     FFTWComplex& operator=(TinyVector<T, 2> const & o)
00311     {
00312         data_[0] = (Real)o[0];
00313         data_[1] = (Real)o[1];
00314         return *this;
00315     }
00316 
00317         /** Assignment.
00318         */
00319     template <class T>
00320     FFTWComplex& operator=(std::complex<T> const & o)
00321     {
00322         data_[0] = (Real)o.real();
00323         data_[1] = (Real)o.imag();
00324         return *this;
00325     }
00326 
00327     reference re()
00328         { return data_[0]; }
00329 
00330     const_reference re() const
00331         { return data_[0]; }
00332 
00333     reference real()
00334         { return data_[0]; }
00335 
00336     const_reference real() const
00337         { return data_[0]; }
00338 
00339     reference im()
00340         { return data_[1]; }
00341 
00342     const_reference im() const
00343         { return data_[1]; }
00344 
00345     reference imag()
00346         { return data_[1]; }
00347 
00348     const_reference imag() const
00349         { return data_[1]; }
00350 
00351         /** Unary negation.
00352         */
00353     FFTWComplex operator-() const
00354         { return FFTWComplex(-data_[0], -data_[1]); }
00355 
00356         /** Squared magnitude x*conj(x)
00357         */
00358     SquaredNormType squaredMagnitude() const
00359         { return data_[0]*data_[0]+data_[1]*data_[1]; }
00360 
00361         /** Magnitude (length of radius vector).
00362         */
00363     NormType magnitude() const
00364         { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
00365 
00366         /** Phase angle.
00367         */
00368     value_type phase() const
00369         { return VIGRA_CSTD::atan2(data_[1], data_[0]); }
00370 
00371         /** Access components as if number were a vector.
00372         */
00373     reference operator[](int i)
00374         { return data_[i]; }
00375 
00376         /** Read components as if number were a vector.
00377         */
00378     const_reference operator[](int i) const
00379         { return data_[i]; }
00380 
00381         /** Length of complex number (always 2).
00382         */
00383     int size() const
00384         { return 2; }
00385 
00386     iterator begin()
00387         { return data_; }
00388 
00389     iterator end()
00390         { return data_ + 2; }
00391 
00392     const_iterator begin() const
00393         { return data_; }
00394 
00395     const_iterator end() const
00396         { return data_ + 2; }
00397 
00398   private:
00399     complex_type data_;
00400 };
00401 
00402 /********************************************************/
00403 /*                                                      */
00404 /*                     FFTWComplexTraits                */
00405 /*                                                      */
00406 /********************************************************/
00407 
00408 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex
00409 
00410     The numeric and promote traits for fftw_complex and FFTWComplex follow
00411     the general specifications for \ref NumericPromotionTraits and
00412     \ref AlgebraicField. They are explicitly specialized for the types
00413     involved:
00414 
00415     \code
00416 
00417     template<>
00418     struct NumericTraits<fftw_complex>
00419     {
00420         typedef fftw_complex Promote;
00421         typedef fftw_complex RealPromote;
00422         typedef fftw_complex ComplexPromote;
00423         typedef fftw_real    ValueType;
00424 
00425         typedef VigraFalseType isIntegral;
00426         typedef VigraFalseType isScalar;
00427         typedef VigraFalseType isOrdered;
00428         typedef VigraTrueType  isComplex;
00429 
00430         // etc.
00431     };
00432 
00433     template<class Real>
00434     struct NumericTraits<FFTWComplex<Real> >
00435     {
00436         typedef FFTWComplex<Real> Promote;
00437         typedef FFTWComplex<Real> RealPromote;
00438         typedef FFTWComplex<Real> ComplexPromote;
00439         typedef Real              ValueType;
00440 
00441         typedef VigraFalseType isIntegral;
00442         typedef VigraFalseType isScalar;
00443         typedef VigraFalseType isOrdered;
00444         typedef VigraTrueType  isComplex;
00445 
00446         // etc.
00447     };
00448 
00449     template<>
00450     struct NormTraits<fftw_complex>
00451     {
00452         typedef fftw_complex Type;
00453         typedef fftw_real    SquaredNormType;
00454         typedef fftw_real    NormType;
00455     };
00456 
00457     template<class Real>
00458     struct NormTraits<FFTWComplex>
00459     {
00460         typedef FFTWComplex Type;
00461         typedef fftw_real   SquaredNormType;
00462         typedef fftw_real   NormType;
00463     };
00464 
00465     template <>
00466     struct PromoteTraits<fftw_complex, fftw_complex>
00467     {
00468         typedef fftw_complex Promote;
00469     };
00470 
00471     template <>
00472     struct PromoteTraits<fftw_complex, double>
00473     {
00474         typedef fftw_complex Promote;
00475     };
00476 
00477     template <>
00478     struct PromoteTraits<double, fftw_complex>
00479     {
00480         typedef fftw_complex Promote;
00481     };
00482 
00483     template <class Real>
00484     struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
00485     {
00486         typedef FFTWComplex<Real> Promote;
00487     };
00488 
00489     template <class Real>
00490     struct PromoteTraits<FFTWComplex<Real>, double>
00491     {
00492         typedef FFTWComplex<Real> Promote;
00493     };
00494 
00495     template <class Real>
00496     struct PromoteTraits<double, FFTWComplex<Real> >
00497     {
00498         typedef FFTWComplex<Real> Promote;
00499     };
00500     \endcode
00501 
00502     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
00503     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
00504     Namespace: vigra
00505 
00506 */
00507 template<>
00508 struct NumericTraits<fftw_complex>
00509 {
00510     typedef fftw_complex Type;
00511     typedef fftw_complex Promote;
00512     typedef fftw_complex RealPromote;
00513     typedef fftw_complex ComplexPromote;
00514     typedef fftw_real    ValueType;
00515 
00516     typedef VigraFalseType isIntegral;
00517     typedef VigraFalseType isScalar;
00518     typedef NumericTraits<fftw_real>::isSigned isSigned;
00519     typedef VigraFalseType isOrdered;
00520     typedef VigraTrueType  isComplex;
00521 
00522     static FFTWComplex<> zero() { return FFTWComplex<>(0.0, 0.0); }
00523     static FFTWComplex<> one() { return FFTWComplex<>(1.0, 0.0); }
00524     static FFTWComplex<> nonZero() { return one(); }
00525 
00526     static const Promote & toPromote(const Type & v) { return v; }
00527     static const RealPromote & toRealPromote(const Type & v) { return v; }
00528     static const Type & fromPromote(const Promote & v) { return v; }
00529     static const Type & fromRealPromote(const RealPromote & v) { return v; }
00530 };
00531 
00532 template<class Real>
00533 struct NumericTraits<FFTWComplex<Real> >
00534 {
00535     typedef FFTWComplex<Real> Type;
00536     typedef FFTWComplex<Real> Promote;
00537     typedef FFTWComplex<Real> RealPromote;
00538     typedef FFTWComplex<Real> ComplexPromote;
00539     typedef typename Type::value_type ValueType;
00540 
00541     typedef VigraFalseType isIntegral;
00542     typedef VigraFalseType isScalar;
00543     typedef typename NumericTraits<ValueType>::isSigned isSigned;
00544     typedef VigraFalseType isOrdered;
00545     typedef VigraTrueType  isComplex;
00546 
00547     static Type zero() { return Type(0.0, 0.0); }
00548     static Type one() { return Type(1.0, 0.0); }
00549     static Type nonZero() { return one(); }
00550 
00551     static const Promote & toPromote(const Type & v) { return v; }
00552     static const RealPromote & toRealPromote(const Type & v) { return v; }
00553     static const Type & fromPromote(const Promote & v) { return v; }
00554     static const Type & fromRealPromote(const RealPromote & v) { return v; }
00555 };
00556 
00557 template<>
00558 struct NormTraits<fftw_complex>
00559 {
00560     typedef fftw_complex Type;
00561     typedef fftw_real    SquaredNormType;
00562     typedef fftw_real    NormType;
00563 };
00564 
00565 template<class Real>
00566 struct NormTraits<FFTWComplex<Real> >
00567 {
00568     typedef FFTWComplex<Real>  Type;
00569     typedef typename Type::SquaredNormType   SquaredNormType;
00570     typedef typename Type::NormType   NormType;
00571 };
00572 
00573 template <>
00574 struct PromoteTraits<fftw_complex, fftw_complex>
00575 {
00576     typedef fftw_complex Promote;
00577 };
00578 
00579 template <>
00580 struct PromoteTraits<fftw_complex, double>
00581 {
00582     typedef fftw_complex Promote;
00583 };
00584 
00585 template <>
00586 struct PromoteTraits<double, fftw_complex>
00587 {
00588     typedef fftw_complex Promote;
00589 };
00590 
00591 template <class Real>
00592 struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
00593 {
00594     typedef FFTWComplex<Real> Promote;
00595 };
00596 
00597 template <class Real>
00598 struct PromoteTraits<FFTWComplex<Real>, double>
00599 {
00600     typedef FFTWComplex<Real> Promote;
00601 };
00602 
00603 template <class Real>
00604 struct PromoteTraits<double, FFTWComplex<Real> >
00605 {
00606     typedef FFTWComplex<Real> Promote;
00607 };
00608 
00609 template<class T>
00610 struct CanSkipInitialization<std::complex<T> >
00611 {
00612     typedef typename CanSkipInitialization<T>::type type;
00613     static const bool value = type::asBool;
00614 };
00615 
00616 template<class Real>
00617 struct CanSkipInitialization<FFTWComplex<Real> >
00618 {
00619     typedef typename CanSkipInitialization<Real>::type type;
00620     static const bool value = type::asBool;
00621 };
00622 
00623 namespace multi_math {
00624 
00625 template <class ARG>
00626 struct MultiMathOperand;
00627 
00628 template <class Real>
00629 struct MultiMathOperand<FFTWComplex<Real> >
00630 {
00631     typedef MultiMathOperand<FFTWComplex<Real> > AllowOverload;
00632     typedef FFTWComplex<Real> result_type;
00633     
00634     static const int ndim = 0;
00635     
00636     MultiMathOperand(FFTWComplex<Real> const & v)
00637     : v_(v)
00638     {}
00639     
00640     template <class SHAPE>
00641     bool checkShape(SHAPE const &) const
00642     {
00643         return true;
00644     }
00645     
00646     template <class SHAPE>
00647     FFTWComplex<Real> const & operator[](SHAPE const &) const
00648     {
00649         return v_;
00650     }
00651 
00652     void inc(unsigned int /*LEVEL*/) const
00653     {}
00654 
00655     void reset(unsigned int /*LEVEL*/) const
00656     {}
00657     
00658     FFTWComplex<Real> const & operator*() const
00659     {
00660         return v_;
00661     }
00662     
00663     FFTWComplex<Real> v_;
00664 };
00665 
00666 } // namespace multi_math
00667 
00668 template<class Ty>
00669 class FFTWAllocator
00670 {
00671   public:
00672     typedef std::size_t size_type;
00673     typedef std::ptrdiff_t difference_type;
00674     typedef Ty *pointer;
00675     typedef const Ty *const_pointer;
00676     typedef Ty& reference;
00677     typedef const Ty& const_reference;
00678     typedef Ty value_type;
00679     
00680     pointer address(reference val) const
00681         { return &val; }
00682         
00683     const_pointer address(const_reference val) const
00684         { return &val; }
00685         
00686     template<class Other>
00687     struct rebind
00688     {
00689         typedef FFTWAllocator<Other> other;
00690     };
00691     
00692     FFTWAllocator() throw()
00693     {}
00694     
00695     template<class Other>
00696     FFTWAllocator(const FFTWAllocator<Other>& right) throw()
00697     {}
00698     
00699     template<class Other>
00700     FFTWAllocator& operator=(const FFTWAllocator<Other>& right)
00701     {
00702         return *this;
00703     }
00704     
00705     pointer allocate(size_type count, void * = 0)
00706     {
00707         return (pointer)fftw_malloc(count * sizeof(Ty));
00708     }
00709     
00710     void deallocate(pointer ptr, size_type count)
00711     {
00712         fftw_free(ptr);
00713     }
00714     
00715     void construct(pointer ptr, const Ty& val)
00716     {
00717         new(ptr) Ty(val);
00718         
00719     }
00720     
00721     void destroy(pointer ptr)
00722     {
00723         ptr->~Ty();
00724     }
00725     
00726     size_type max_size() const throw()
00727     {
00728         return NumericTraits<std::ptrdiff_t>::max() / sizeof(Ty);
00729     }
00730 };
00731 
00732 } // namespace vigra
00733 
00734 namespace std {
00735 
00736 template<class Real>
00737 class allocator<vigra::FFTWComplex<Real> >
00738 {
00739   public:
00740     typedef vigra::FFTWComplex<Real> value_type;
00741     typedef size_t size_type;
00742     typedef std::ptrdiff_t difference_type;
00743     typedef value_type *pointer;
00744     typedef const value_type *const_pointer;
00745     typedef value_type& reference;
00746     typedef const value_type& const_reference;
00747 
00748     pointer address(reference val) const
00749         { return &val; }
00750         
00751     const_pointer address(const_reference val) const
00752         { return &val; }
00753         
00754     template<class Other>
00755     struct rebind
00756     {
00757         typedef allocator<Other> other;
00758     };
00759     
00760     allocator() throw()
00761     {}
00762     
00763     template<class Other>
00764     allocator(const allocator<Other>& right) throw()
00765     {}
00766     
00767     template<class Other>
00768     allocator& operator=(const allocator<Other>& right)
00769     {
00770         return *this;
00771     }
00772     
00773     pointer allocate(size_type count, void * = 0)
00774     {
00775         return (pointer)fftw_malloc(count * sizeof(value_type));
00776     }
00777 
00778     void deallocate(pointer ptr, size_type /*count*/)
00779     {
00780         fftw_free(ptr);
00781     }
00782     
00783     void construct(pointer ptr, const value_type& val)
00784     {
00785         new(ptr) value_type(val);
00786         
00787     }
00788     
00789     void destroy(pointer ptr)
00790     {
00791         ptr->~value_type();
00792     }
00793     
00794     size_type max_size() const throw()
00795     {
00796         return vigra::NumericTraits<std::ptrdiff_t>::max() / sizeof(value_type);
00797     }
00798 };
00799 
00800 } // namespace std
00801 
00802 namespace vigra {
00803 
00804 /********************************************************/
00805 /*                                                      */
00806 /*                    FFTWComplex Operations            */
00807 /*                                                      */
00808 /********************************************************/
00809 
00810 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex
00811 
00812     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
00813     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
00814 
00815     These functions fulfill the requirements of an Algebraic Field.
00816     Return types are determined according to \ref FFTWComplexTraits.
00817 
00818     Namespace: vigra
00819     <p>
00820 
00821  */
00822 //@{
00823     /// equal
00824 template <class R>
00825 inline bool operator ==(FFTWComplex<R> const &a, const FFTWComplex<R> &b) {
00826     return a.re() == b.re() && a.im() == b.im();
00827 }
00828 
00829 template <class R>
00830 inline bool operator ==(FFTWComplex<R> const &a,  double b) {
00831     return a.re() == b && a.im() == 0.0;
00832 }
00833 
00834 template <class R>
00835 inline bool operator ==(double a, const FFTWComplex<R> &b) {
00836     return a == b.re() && 0.0 == b.im();
00837 }
00838 
00839     /// not equal
00840 template <class R>
00841 inline bool operator !=(FFTWComplex<R> const &a, const FFTWComplex<R> &b) {
00842     return a.re() != b.re() || a.im() != b.im();
00843 }
00844 
00845     /// not equal
00846 template <class R>
00847 inline bool operator !=(FFTWComplex<R> const &a,  double b) {
00848     return a.re() != b || a.im() != 0.0;
00849 }
00850 
00851     /// not equal
00852 template <class R>
00853 inline bool operator !=(double a, const FFTWComplex<R> &b) {
00854     return a != b.re() || 0.0 != b.im();
00855 }
00856 
00857     /// add-assignment
00858 template <class R>
00859 inline FFTWComplex<R> & operator +=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
00860     a.re() += b.re();
00861     a.im() += b.im();
00862     return a;
00863 }
00864 
00865     /// subtract-assignment
00866 template <class R>
00867 inline FFTWComplex<R> & operator -=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
00868     a.re() -= b.re();
00869     a.im() -= b.im();
00870     return a;
00871 }
00872 
00873     /// multiply-assignment
00874 template <class R>
00875 inline FFTWComplex<R> & operator *=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
00876     typename FFTWComplex<R>::value_type t = a.re()*b.re()-a.im()*b.im();
00877     a.im() = a.re()*b.im()+a.im()*b.re();
00878     a.re() = t;
00879     return a;
00880 }
00881 
00882     /// divide-assignment
00883 template <class R>
00884 inline FFTWComplex<R> & operator /=(FFTWComplex<R> &a, const FFTWComplex<R> &b) {
00885     typename FFTWComplex<R>::value_type sm = b.squaredMagnitude();
00886     typename FFTWComplex<R>::value_type t = (a.re()*b.re()+a.im()*b.im())/sm;
00887     a.im() = (b.re()*a.im()-a.re()*b.im())/sm;
00888     a.re() = t;
00889     return a;
00890 }
00891 
00892     /// add-assignment with scalar double
00893 template <class R>
00894 inline FFTWComplex<R> & operator +=(FFTWComplex<R> &a, double b) {
00895     a.re() += (R)b;
00896     return a;
00897 }
00898 
00899     /// subtract-assignment with scalar double
00900 template <class R>
00901 inline FFTWComplex<R> & operator -=(FFTWComplex<R> &a, double b) {
00902     a.re() -= (R)b;
00903     return a;
00904 }
00905 
00906     /// multiply-assignment with scalar double
00907 template <class R>
00908 inline FFTWComplex<R> & operator *=(FFTWComplex<R> &a, double b) {
00909     a.re() *= (R)b;
00910     a.im() *= (R)b;
00911     return a;
00912 }
00913 
00914     /// divide-assignment with scalar double
00915 template <class R>
00916 inline FFTWComplex<R> & operator /=(FFTWComplex<R> &a, double b) {
00917     a.re() /= (R)b;
00918     a.im() /= (R)b;
00919     return a;
00920 }
00921 
00922     /// addition
00923 template <class R>
00924 inline FFTWComplex<R> operator +(FFTWComplex<R> a, const FFTWComplex<R> &b) {
00925     a += b;
00926     return a;
00927 }
00928 
00929     /// right addition with scalar double
00930 template <class R>
00931 inline FFTWComplex<R> operator +(FFTWComplex<R> a, double b) {
00932     a += b;
00933     return a;
00934 }
00935 
00936     /// left addition with scalar double
00937 template <class R>
00938 inline FFTWComplex<R> operator +(double a, FFTWComplex<R> b) {
00939     b += a;
00940     return b;
00941 }
00942 
00943     /// subtraction
00944 template <class R>
00945 inline FFTWComplex<R> operator -(FFTWComplex<R> a, const FFTWComplex<R> &b) {
00946     a -= b;
00947     return a;
00948 }
00949 
00950     /// right subtraction with scalar double
00951 template <class R>
00952 inline FFTWComplex<R> operator -(FFTWComplex<R> a, double b) {
00953     a -= b;
00954     return a;
00955 }
00956 
00957     /// left subtraction with scalar double
00958 template <class R>
00959 inline FFTWComplex<R> operator -(double a, FFTWComplex<R> const & b) {
00960     return (-b) += a;
00961 }
00962 
00963     /// multiplication
00964 template <class R>
00965 inline FFTWComplex<R> operator *(FFTWComplex<R> a, const FFTWComplex<R> &b) {
00966     a *= b;
00967     return a;
00968 }
00969 
00970     /// right multiplication with scalar double
00971 template <class R>
00972 inline FFTWComplex<R> operator *(FFTWComplex<R> a, double b) {
00973     a *= b;
00974     return a;
00975 }
00976 
00977     /// left multiplication with scalar double
00978 template <class R>
00979 inline FFTWComplex<R> operator *(double a, FFTWComplex<R> b) {
00980     b *= a;
00981     return b;
00982 }
00983 
00984     /// division
00985 template <class R>
00986 inline FFTWComplex<R> operator /(FFTWComplex<R> a, const FFTWComplex<R> &b) {
00987     a /= b;
00988     return a;
00989 }
00990 
00991     /// right division with scalar double
00992 template <class R>
00993 inline FFTWComplex<R> operator /(FFTWComplex<R> a, double b) {
00994     a /= b;
00995     return a;
00996 }
00997 
00998 using VIGRA_CSTD::abs;
00999 
01000     /// absolute value (= magnitude)
01001 template <class R>
01002 inline typename FFTWComplex<R>::NormType abs(const FFTWComplex<R> &a)
01003 {
01004     return a.magnitude();
01005 }
01006 
01007     /// pahse
01008 template <class R>
01009 inline R arg(const FFTWComplex<R> &a)
01010 {
01011     return a.phase();
01012 }
01013 
01014     /// real part
01015 template <class R>
01016 inline R real(const FFTWComplex<R> &a)
01017 {
01018     return a.real();
01019 }
01020 
01021     /// imaginary part
01022 template <class R>
01023 inline R imag(const FFTWComplex<R> &a)
01024 {
01025     return a.imag();
01026 }
01027 
01028     /// complex conjugate
01029 template <class R>
01030 inline FFTWComplex<R> conj(const FFTWComplex<R> &a)
01031 {
01032     return FFTWComplex<R>(a.re(), -a.im());
01033 }
01034 
01035     /// norm (= magnitude)
01036 template <class R>
01037 inline typename FFTWComplex<R>::NormType norm(const FFTWComplex<R> &a)
01038 {
01039     return a.magnitude();
01040 }
01041 
01042     /// squared norm (= squared magnitude)
01043 template <class R>
01044 inline typename FFTWComplex<R>::SquaredNormType squaredNorm(const FFTWComplex<R> &a)
01045 {
01046     return a.squaredMagnitude();
01047 }
01048 
01049 #define VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(fct) \
01050 template <class R> \
01051 inline FFTWComplex<R> fct(const FFTWComplex<R> &a) \
01052 { \
01053     return std::fct(reinterpret_cast<std::complex<R> const &>(a)); \
01054 }
01055 
01056 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cos)
01057 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cosh)
01058 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(exp)
01059 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(log)
01060 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(log10)
01061 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sin)
01062 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sinh)
01063 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sqrt)
01064 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tan)
01065 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tanh)
01066 
01067 #undef VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION
01068 
01069 template <class R>
01070 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, int e)
01071 {
01072     return std::pow(reinterpret_cast<std::complex<R> const &>(a), e);
01073 }
01074 
01075 template <class R>
01076 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, R const & e)
01077 {
01078     return std::pow(reinterpret_cast<std::complex<R> const &>(a), e);
01079 }
01080 
01081 template <class R>
01082 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, const FFTWComplex<R> & e)
01083 {
01084     return std::pow(reinterpret_cast<std::complex<R> const &>(a), 
01085                      reinterpret_cast<std::complex<R> const &>(e));
01086 }
01087 
01088 template <class R>
01089 inline FFTWComplex<R> pow(R const & a, const FFTWComplex<R> &e)
01090 {
01091     return std::pow(a, reinterpret_cast<std::complex<R> const &>(e));
01092 }
01093 
01094 //@}
01095 
01096 } // namespace vigra
01097 
01098 namespace std {
01099 
01100 template <class Real>
01101 ostream & operator<<(ostream & s, vigra::FFTWComplex<Real> const & v)
01102 {
01103     s << std::complex<Real>(v.re(), v.im());
01104     return s;
01105 }
01106 
01107 } // namespace std
01108 
01109 namespace vigra {
01110 
01111 /** \addtogroup StandardImageTypes
01112 */
01113 //@{
01114 
01115 /********************************************************/
01116 /*                                                      */
01117 /*                      FFTWRealImage                   */
01118 /*                                                      */
01119 /********************************************************/
01120 
01121     /** Float (<tt>fftw_real</tt>) image.
01122 
01123         The type <tt>fftw_real</tt> is defined as <tt>double</tt> (in FFTW 2 it used to be
01124         either <tt>float</tt> or <tt>double</tt>, as specified during compilation of FFTW).
01125         FFTWRealImage uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
01126         their const counterparts to access the data.
01127 
01128         <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01129         <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01130         Namespace: vigra
01131     */
01132 typedef BasicImage<fftw_real> FFTWRealImage;
01133 
01134 /********************************************************/
01135 /*                                                      */
01136 /*                     FFTWComplexImage                 */
01137 /*                                                      */
01138 /********************************************************/
01139 
01140 template<class R>
01141 struct IteratorTraits<
01142         BasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> >
01143 {
01144     typedef FFTWComplex<R> Type;
01145     typedef BasicImageIterator<Type, Type **>       Iterator;
01146     typedef Iterator                                iterator;
01147     typedef BasicImageIterator<Type, Type **>       mutable_iterator;
01148     typedef ConstBasicImageIterator<Type, Type **>  const_iterator;
01149     typedef typename iterator::iterator_category    iterator_category;
01150     typedef typename iterator::value_type           value_type;
01151     typedef typename iterator::reference            reference;
01152     typedef typename iterator::index_reference      index_reference;
01153     typedef typename iterator::pointer              pointer;
01154     typedef typename iterator::difference_type      difference_type;
01155     typedef typename iterator::row_iterator         row_iterator;
01156     typedef typename iterator::column_iterator      column_iterator;
01157     typedef VectorAccessor<Type>                    default_accessor;
01158     typedef VectorAccessor<Type>                    DefaultAccessor;
01159     typedef VigraTrueType                           hasConstantStrides;
01160 };
01161 
01162 template<class R>
01163 struct IteratorTraits<
01164         ConstBasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> >
01165 {
01166     typedef FFTWComplex<R> Type;
01167     typedef ConstBasicImageIterator<Type, Type **> Iterator;
01168     typedef Iterator                               iterator;
01169     typedef BasicImageIterator<Type, Type **>      mutable_iterator;
01170     typedef ConstBasicImageIterator<Type, Type **> const_iterator;
01171     typedef typename iterator::iterator_category   iterator_category;
01172     typedef typename iterator::value_type          value_type;
01173     typedef typename iterator::reference           reference;
01174     typedef typename iterator::index_reference     index_reference;
01175     typedef typename iterator::pointer             pointer;
01176     typedef typename iterator::difference_type     difference_type;
01177     typedef typename iterator::row_iterator        row_iterator;
01178     typedef typename iterator::column_iterator     column_iterator;
01179     typedef VectorAccessor<Type>                   default_accessor;
01180     typedef VectorAccessor<Type>                   DefaultAccessor;
01181     typedef VigraTrueType                          hasConstantStrides;
01182 };
01183 
01184     /** Complex (FFTWComplex) image.
01185         It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
01186         their const counterparts to access the data.
01187 
01188         <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01189         <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01190         Namespace: vigra
01191     */
01192 typedef BasicImage<FFTWComplex<> > FFTWComplexImage;
01193 
01194 //@}
01195 
01196 /********************************************************/
01197 /*                                                      */
01198 /*                  FFTWComplex-Accessors               */
01199 /*                                                      */
01200 /********************************************************/
01201 
01202 /** \addtogroup DataAccessors
01203 */
01204 //@{
01205 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex
01206 
01207     Encapsulate access to pixels of type FFTWComplex
01208 */
01209 //@{
01210     /** Encapsulate access to the the real part of a complex number.
01211 
01212     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01213     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01214     Namespace: vigra
01215     */
01216 template <class Real = double>
01217 class FFTWRealAccessor
01218 {
01219   public:
01220 
01221         /// The accessor's value type.
01222     typedef Real value_type;
01223 
01224         /// Read real part at iterator position.
01225     template <class ITERATOR>
01226     value_type operator()(ITERATOR const & i) const {
01227         return (*i).re();
01228     }
01229 
01230         /// Read real part at offset from iterator position.
01231     template <class ITERATOR, class DIFFERENCE>
01232     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
01233         return i[d].re();
01234     }
01235 
01236         /// Write real part at iterator position from a scalar.
01237     template <class ITERATOR>
01238     void set(value_type const & v, ITERATOR const & i) const {
01239         (*i).re()= v;
01240     }
01241 
01242         /// Write real part at offset from iterator position from a scalar.
01243     template <class ITERATOR, class DIFFERENCE>
01244     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
01245         i[d].re()= v;
01246     }
01247 
01248         /// Write real part at iterator position into a scalar.
01249     template <class R, class ITERATOR>
01250     void set(FFTWComplex<R> const & v, ITERATOR const & i) const {
01251         *i = v.re();
01252     }
01253 
01254         /// Write real part at offset from iterator position into a scalar.
01255     template <class R, class ITERATOR, class DIFFERENCE>
01256     void set(FFTWComplex<R> const & v, ITERATOR const & i, DIFFERENCE d) const {
01257         i[d] = v.re();
01258     }
01259 };
01260 
01261     /** Encapsulate access to the the imaginary part of a complex number.
01262 
01263     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01264     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01265     Namespace: vigra
01266     */
01267 template <class Real = double>
01268 class FFTWImaginaryAccessor
01269 {
01270   public:
01271         /// The accessor's value type.
01272     typedef Real value_type;
01273 
01274         /// Read imaginary part at iterator position.
01275     template <class ITERATOR>
01276     value_type operator()(ITERATOR const & i) const {
01277         return (*i).im();
01278     }
01279 
01280         /// Read imaginary part at offset from iterator position.
01281     template <class ITERATOR, class DIFFERENCE>
01282     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
01283         return i[d].im();
01284     }
01285 
01286         /// Write imaginary part at iterator position from a scalar.
01287     template <class ITERATOR>
01288     void set(value_type const & v, ITERATOR const & i) const {
01289         (*i).im()= v;
01290     }
01291 
01292         /// Write imaginary part at offset from iterator position from a scalar.
01293     template <class ITERATOR, class DIFFERENCE>
01294     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
01295         i[d].im()= v;
01296     }
01297 
01298         /// Write imaginary part at iterator position into a scalar.
01299     template <class R, class ITERATOR>
01300     void set(FFTWComplex<R> const & v, ITERATOR const & i) const {
01301         *i = v.im();
01302     }
01303 
01304         /// Write imaginary part at offset from iterator position into a scalar.
01305     template <class R, class ITERATOR, class DIFFERENCE>
01306     void set(FFTWComplex<R> const & v, ITERATOR const & i, DIFFERENCE d) const {
01307         i[d] = v.im();
01308     }
01309 };
01310 
01311     /** Write a real number into a complex one. The imaginary part is set
01312         to 0.
01313 
01314     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01315     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01316     Namespace: vigra
01317     */
01318 template <class Real = double>
01319 class FFTWWriteRealAccessor
01320 : public FFTWRealAccessor<Real>
01321 {
01322   public:
01323         /// The accessor's value type.
01324     typedef Real value_type;
01325 
01326         /** Write real number at iterator position. Set imaginary part
01327             to 0.
01328         */
01329     template <class ITERATOR>
01330     void set(value_type const & v, ITERATOR const & i) const {
01331         (*i).re()= v;
01332         (*i).im()= 0;
01333     }
01334 
01335         /** Write real number at offset from iterator position. Set imaginary part
01336             to 0.
01337         */
01338     template <class ITERATOR, class DIFFERENCE>
01339     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
01340         i[d].re()= v;
01341         i[d].im()= 0;
01342     }
01343 };
01344 
01345     /** Calculate squared magnitude of complex number on the fly.
01346 
01347     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01348     Namespace: vigra
01349     */
01350 template <class Real = double>
01351 class FFTWSquaredMagnitudeAccessor
01352 {
01353   public:
01354         /// The accessor's value type.
01355     typedef Real value_type;
01356 
01357         /// Read squared magnitude at iterator position.
01358     template <class ITERATOR>
01359     value_type operator()(ITERATOR const & i) const {
01360         return (*i).squaredMagnitude();
01361     }
01362 
01363         /// Read squared magnitude at offset from iterator position.
01364     template <class ITERATOR, class DIFFERENCE>
01365     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
01366         return (i[d]).squaredMagnitude();
01367     }
01368 };
01369 
01370     /** Calculate magnitude of complex number on the fly.
01371 
01372     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01373     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01374     Namespace: vigra
01375     */
01376 template <class Real = double>
01377 class FFTWMagnitudeAccessor
01378 {
01379   public:
01380         /// The accessor's value type.
01381     typedef Real value_type;
01382 
01383         /// Read magnitude at iterator position.
01384     template <class ITERATOR>
01385     value_type operator()(ITERATOR const & i) const {
01386         return (*i).magnitude();
01387     }
01388 
01389         /// Read magnitude at offset from iterator position.
01390     template <class ITERATOR, class DIFFERENCE>
01391     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
01392         return (i[d]).magnitude();
01393     }
01394 };
01395 
01396     /** Calculate phase of complex number on the fly.
01397 
01398     <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br>
01399     <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br>
01400     Namespace: vigra
01401     */
01402 template <class Real = double>
01403 class FFTWPhaseAccessor
01404 {
01405   public:
01406         /// The accessor's value type.
01407     typedef Real value_type;
01408 
01409         /// Read phase at iterator position.
01410     template <class ITERATOR>
01411     value_type operator()(ITERATOR const & i) const {
01412         return (*i).phase();
01413     }
01414 
01415         /// Read phase at offset from iterator position.
01416     template <class ITERATOR, class DIFFERENCE>
01417     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
01418         return (i[d]).phase();
01419     }
01420 };
01421 
01422 //@}
01423 //@}
01424 
01425 /********************************************************/
01426 /*                                                      */
01427 /*                    Fourier Transform                 */
01428 /*                                                      */
01429 /********************************************************/
01430 
01431 /* \addtogroup FourierTransform Fast Fourier Transform
01432 
01433     This documentation describes the VIGRA interface to FFTW version 3. The interface
01434     to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated.
01435 
01436     VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
01437     Transform</a> package to perform Fourier transformations. VIGRA
01438     provides a wrapper for FFTW's complex number type (FFTWComplex),
01439     but FFTW's functions are used verbatim. If the image is stored as
01440     a FFTWComplexImage, the simplest call to an FFT function is like this:
01441 
01442     \code
01443     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
01444     ... // fill image with data
01445 
01446     // create a plan with estimated performance optimization
01447     fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
01448                                 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
01449                                 FFTW_FORWARD, FFTW_ESTIMATE );
01450     // calculate FFT (this can be repeated as often as needed,
01451     //                with fresh data written into the source array)
01452     fftw_execute(forwardPlan);
01453 
01454     // release the plan memory
01455     fftw_destroy_plan(forwardPlan);
01456 
01457     // likewise for the inverse transform
01458     fftw_plan backwardPlan = fftw_plan_dft_2d(height, width,
01459                                  (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(),
01460                                  FFTW_BACKWARD, FFTW_ESTIMATE);
01461     fftw_execute(backwardPlan);
01462     fftw_destroy_plan(backwardPlan);
01463 
01464     // do not forget to normalize the result according to the image size
01465     transformImage(srcImageRange(spatial), destImage(spatial),
01466                    std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height));
01467     \endcode
01468 
01469     Note that in the creation of a plan, the height must be given
01470     first. Note also that <TT>spatial.begin()</TT> may only be passed
01471     to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the
01472     entire image. When you want to restrict operation to an ROI, you
01473     can create a copy of the ROI in an image of appropriate size, or
01474     you may use the Guru interface to FFTW.
01475 
01476     More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
01477 
01478     FFTW produces fourier images that have the DC component (the
01479     origin of the Fourier space) in the upper left corner. Often, one
01480     wants the origin in the center of the image, so that frequencies
01481     always increase towards the border of the image. This can be
01482     achieved by calling \ref moveDCToCenter(). The inverse
01483     transformation is done by \ref moveDCToUpperLeft().
01484 
01485     <b>\#include</b> <vigra/fftw3.hxx><br>
01486     Namespace: vigra
01487 */
01488 
01489 /** \addtogroup FourierTransform Fast Fourier Transform
01490 
01491     VIGRA provides a powerful C++ API for the popular <a href="http://www.fftw.org/">FFTW library</a>
01492     for fast Fourier transforms. There are two versions of the API: an older one based on image 
01493     iterators (and therefore restricted to 2D) and a new one based on \ref MultiArrayView that
01494     works for arbitrary dimensions. In addition, the functions \ref convolveFFT() and 
01495     \ref applyFourierFilter() provide an easy-to-use interface for FFT-based convolution,
01496     a major application of Fourier transforms.
01497 */
01498 //@{
01499 
01500 /********************************************************/
01501 /*                                                      */
01502 /*                     moveDCToCenter                   */
01503 /*                                                      */
01504 /********************************************************/
01505 
01506 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
01507           in the image center.
01508 
01509     FFTW produces fourier images where the DC component (origin of
01510     fourier space) is located in the upper left corner of the
01511     image. The quadrants are placed like this (using a 4x4 image for
01512     example):
01513 
01514     \code
01515             DC 4 3 3
01516              4 4 3 3
01517              1 1 2 2
01518              1 1 2 2
01519     \endcode
01520 
01521     After applying the function, the quadrants are at their usual places:
01522 
01523     \code
01524             2 2  1 1
01525             2 2  1 1
01526             3 3 DC 4
01527             3 3  4 4
01528     \endcode
01529 
01530     This transformation can be reversed by \ref moveDCToUpperLeft().
01531     Note that the 2D versions of this transformation must not be executed in place - input
01532     and output images must be different. In contrast, the nD version (with MultiArrayView
01533     argument) always works in-place.
01534 
01535     <b> Declarations:</b>
01536 
01537     use MultiArrayView (this works in-place, with arbitrary dimension N):
01538     \code
01539     namespace vigra {
01540         template <unsigned int N, class T, class Stride>
01541         void moveDCToCenter(MultiArrayView<N, T, Stride> a);
01542     }
01543     \endcode
01544 
01545     pass iterators explicitly (2D only, not in-place):
01546     \code
01547     namespace vigra {
01548         template <class SrcImageIterator, class SrcAccessor,
01549                   class DestImageIterator, class DestAccessor>
01550         void moveDCToCenter(SrcImageIterator src_upperleft,
01551                                SrcImageIterator src_lowerright, SrcAccessor sa,
01552                                DestImageIterator dest_upperleft, DestAccessor da);
01553     }
01554     \endcode
01555 
01556 
01557     use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place):
01558     \code
01559     namespace vigra {
01560         template <class SrcImageIterator, class SrcAccessor,
01561                   class DestImageIterator, class DestAccessor>
01562         void moveDCToCenter(
01563             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01564             pair<DestImageIterator, DestAccessor> dest);
01565     }
01566     \endcode
01567 
01568     <b> Usage:</b>
01569 
01570     <b>\#include</b> <vigra/fftw3.hxx> (for 2D variants) <br>
01571     <b>\#include</b> <vigra/multi_fft.hxx> (for nD variants) <br>
01572     Namespace: vigra
01573 
01574     \code
01575     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
01576     ... // fill image with data
01577 
01578     // create a plan with estimated performance optimization
01579     fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
01580                                 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
01581                                 FFTW_FORWARD, FFTW_ESTIMATE );
01582     // calculate FFT
01583     fftw_execute(forwardPlan);
01584 
01585     vigra::FFTWComplexImage rearrangedFourier(width, height);
01586     moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
01587 
01588     // delete the plan
01589     fftw_destroy_plan(forwardPlan);
01590     \endcode
01591 */
01592 doxygen_overloaded_function(template <...> void moveDCToCenter)
01593 
01594 template <class SrcImageIterator, class SrcAccessor,
01595           class DestImageIterator, class DestAccessor>
01596 void moveDCToCenter(SrcImageIterator src_upperleft,
01597                                SrcImageIterator src_lowerright, SrcAccessor sa,
01598                                DestImageIterator dest_upperleft, DestAccessor da)
01599 {
01600     int w = int(src_lowerright.x - src_upperleft.x);
01601     int h = int(src_lowerright.y - src_upperleft.y);
01602     int w1 = w/2;
01603     int h1 = h/2;
01604     int w2 = (w+1)/2;
01605     int h2 = (h+1)/2;
01606 
01607     // 2. Quadrant  zum 4.
01608     copyImage(srcIterRange(src_upperleft,
01609                            src_upperleft  + Diff2D(w2, h2), sa),
01610               destIter    (dest_upperleft + Diff2D(w1, h1), da));
01611 
01612     // 4. Quadrant zum 2.
01613     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
01614                            src_lowerright, sa),
01615               destIter    (dest_upperleft, da));
01616 
01617     // 1. Quadrant zum 3.
01618     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
01619                            src_upperleft  + Diff2D(w,  h2), sa),
01620               destIter    (dest_upperleft + Diff2D(0,  h1), da));
01621 
01622     // 3. Quadrant zum 1.
01623     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
01624                            src_upperleft  + Diff2D(w2, h), sa),
01625               destIter    (dest_upperleft + Diff2D(w1, 0), da));
01626 }
01627 
01628 template <class SrcImageIterator, class SrcAccessor,
01629           class DestImageIterator, class DestAccessor>
01630 inline void moveDCToCenter(
01631     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01632     pair<DestImageIterator, DestAccessor> dest)
01633 {
01634     moveDCToCenter(src.first, src.second, src.third,
01635                           dest.first, dest.second);
01636 }
01637 
01638 /********************************************************/
01639 /*                                                      */
01640 /*                   moveDCToUpperLeft                  */
01641 /*                                                      */
01642 /********************************************************/
01643 
01644 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
01645           in the image's upper left.
01646 
01647      This function is the inversion of \ref moveDCToCenter(). See there
01648      for a detailed description and usage examples.
01649 
01650      <b> Declarations:</b>
01651 
01652     use MultiArrayView (this works in-place, with arbitrary dimension N):
01653     \code
01654     namespace vigra {
01655         template <unsigned int N, class T, class Stride>
01656         void moveDCToUpperLeft(MultiArrayView<N, T, Stride> a);
01657     }
01658     \endcode
01659 
01660     pass iterators explicitly (2D only, not in-place):
01661      \code
01662         namespace vigra {
01663             template <class SrcImageIterator, class SrcAccessor,
01664                       class DestImageIterator, class DestAccessor>
01665             void moveDCToUpperLeft(SrcImageIterator src_upperleft,
01666                                    SrcImageIterator src_lowerright, SrcAccessor sa,
01667                                    DestImageIterator dest_upperleft, DestAccessor da);
01668         }
01669      \endcode
01670 
01671 
01672      use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place):
01673      \code
01674         namespace vigra {
01675             template <class SrcImageIterator, class SrcAccessor,
01676                       class DestImageIterator, class DestAccessor>
01677             void moveDCToUpperLeft(
01678                 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01679                 pair<DestImageIterator, DestAccessor> dest);
01680         }
01681      \endcode
01682 */
01683 doxygen_overloaded_function(template <...> void moveDCToUpperLeft)
01684 
01685 template <class SrcImageIterator, class SrcAccessor,
01686           class DestImageIterator, class DestAccessor>
01687 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
01688                                SrcImageIterator src_lowerright, SrcAccessor sa,
01689                                DestImageIterator dest_upperleft, DestAccessor da)
01690 {
01691     int w = int(src_lowerright.x - src_upperleft.x);
01692     int h = int(src_lowerright.y - src_upperleft.y);
01693     int w2 = w/2;
01694     int h2 = h/2;
01695     int w1 = (w+1)/2;
01696     int h1 = (h+1)/2;
01697 
01698     // 2. Quadrant  zum 4.
01699     copyImage(srcIterRange(src_upperleft,
01700                            src_upperleft  + Diff2D(w2, h2), sa),
01701               destIter    (dest_upperleft + Diff2D(w1, h1), da));
01702 
01703     // 4. Quadrant zum 2.
01704     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
01705                            src_lowerright, sa),
01706               destIter    (dest_upperleft, da));
01707 
01708     // 1. Quadrant zum 3.
01709     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
01710                            src_upperleft  + Diff2D(w,  h2), sa),
01711               destIter    (dest_upperleft + Diff2D(0,  h1), da));
01712 
01713     // 3. Quadrant zum 1.
01714     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
01715                            src_upperleft  + Diff2D(w2, h), sa),
01716               destIter    (dest_upperleft + Diff2D(w1, 0), da));
01717 }
01718 
01719 template <class SrcImageIterator, class SrcAccessor,
01720           class DestImageIterator, class DestAccessor>
01721 inline void moveDCToUpperLeft(
01722     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01723     pair<DestImageIterator, DestAccessor> dest)
01724 {
01725     moveDCToUpperLeft(src.first, src.second, src.third,
01726                                           dest.first, dest.second);
01727 }
01728 
01729 namespace detail {
01730 
01731 template <class T>
01732 void
01733 fourierTransformImpl(FFTWComplexImage::const_traverser sul,
01734                      FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01735                      FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest, T sign)
01736 {
01737     int w = int(slr.x - sul.x);
01738     int h = int(slr.y - sul.y);
01739 
01740     FFTWComplexImage sworkImage, dworkImage;
01741 
01742     fftw_complex * srcPtr = (fftw_complex *)(&*sul);
01743     fftw_complex * destPtr = (fftw_complex *)(&*dul);
01744 
01745     // test for right memory layout (fftw expects a 2*width*height floats array)
01746     if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1))))
01747     {
01748         sworkImage.resize(w, h);
01749         copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
01750         srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft()));
01751     }
01752     if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
01753     {
01754         dworkImage.resize(w, h);
01755         destPtr = (fftw_complex *)(&(*dworkImage.upperLeft()));
01756     }
01757 
01758     fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
01759     fftw_execute(plan);
01760     fftw_destroy_plan(plan);
01761 
01762     if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
01763     {
01764         copyImage(srcImageRange(dworkImage), destIter(dul, dest));
01765     }
01766 }
01767 
01768 } // namespace detail
01769 
01770 /********************************************************/
01771 /*                                                      */
01772 /*                   fourierTransform                   */
01773 /*                                                      */
01774 /********************************************************/
01775 
01776 /** \brief Compute forward and inverse Fourier transforms.
01777 
01778     The array referring to the spatial domain (i.e. the input in a forward transform, 
01779     and the output in an inverse transform) may be scalar or complex. The array representing
01780     the frequency domain (i.e. output for forward transform, input for inverse transform) 
01781     must always be complex.
01782     
01783     The new implementations (those using MultiArrayView arguments) perform a normalized transform, 
01784     whereas the old ones (using 2D iterators or argument objects) perform an un-normalized 
01785     transform (i.e. the result of the inverse transform is scaled by the number of pixels).
01786 
01787     In general, input and output arrays must have the same shape, with the exception of the 
01788     special <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C 
01789     and C2R modes</a> defined by FFTW.
01790     
01791     The R2C transform reduces the redundancy in the Fourier representation of a real-valued signal:
01792     Since the Fourier representation of a real signal is symmetric, about half of the Fourier coefficients 
01793     can simply be dropped. By convention, this reduction is applied to the first (innermost) dimension, 
01794     such that <tt>fourier.shape(0) == spatial.shape(0)/2 + 1</tt> holds. The correct frequency domain
01795     shape can be conveniently computed by means of the function \ref fftwCorrespondingShapeR2C().
01796     
01797     Note that your program must always link against <tt>libfftw3</tt>. If you want to compute Fourier 
01798     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>).
01799     
01800     The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a>
01801     which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e.
01802     optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning,
01803     you can use the class \ref FFTWPlan.
01804     
01805     <b> Declarations:</b>
01806 
01807     use complex-valued MultiArrayView arguments (this works for arbitrary dimension N):
01808     \code
01809     namespace vigra {
01810         template <unsigned int N, class Real, class C1, class C2>
01811         void 
01812         fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
01813                          MultiArrayView<N, FFTWComplex<Real>, C2> out);
01814 
01815         template <unsigned int N, class Real, class C1, class C2>
01816         void 
01817         fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
01818                                 MultiArrayView<N, FFTWComplex<Real>, C2> out);
01819     }
01820     \endcode
01821 
01822     use real-valued MultiArrayView in the spatial domain, complex-valued MultiArrayView 
01823     in the frequency domain (this works for arbitrary dimension N, and also supports
01824     the R2C and C2R transform, depending on the array shape in the frequency domain):
01825     \code
01826     namespace vigra {
01827         template <unsigned int N, class Real, class C1, class C2>
01828         void 
01829         fourierTransform(MultiArrayView<N, Real, C1> in, 
01830                          MultiArrayView<N, FFTWComplex<Real>, C2> out);
01831 
01832         template <unsigned int N, class Real, class C1, class C2>
01833         void 
01834         fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 
01835                                 MultiArrayView<N, Real, C2> out);
01836     }
01837     \endcode
01838 
01839     pass iterators explicitly (2D only, double only):
01840     \code
01841     namespace vigra {
01842         template <class SrcImageIterator, class SrcAccessor>
01843         void fourierTransform(SrcImageIterator srcUpperLeft,
01844                               SrcImageIterator srcLowerRight, SrcAccessor src,
01845                               FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest);
01846 
01847         void
01848         fourierTransformInverse(FFTWComplexImage::const_traverser sul,
01849                                 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01850                                 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
01851     }
01852     \endcode
01853 
01854     use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, double only):
01855     \code
01856     namespace vigra {
01857         template <class SrcImageIterator, class SrcAccessor>
01858         void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01859                               pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
01860 
01861         void
01862         fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
01863                                        FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
01864                                 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
01865     }
01866     \endcode
01867 
01868     <b> Usage:</b>
01869 
01870     <b>\#include</b> <vigra/fftw3.hxx> (old-style 2D variants)<br>
01871     <b>\#include</b> <vigra/multi_fft.hxx> (new-style nD variants)<br>
01872     Namespace: vigra
01873 
01874     old-style example (using FFTWComplexImage):
01875     \code
01876     // compute complex Fourier transform of a real image, old-style implementation
01877     vigra::DImage src(w, h);
01878     vigra::FFTWComplexImage fourier(w, h);
01879 
01880     fourierTransform(srcImageRange(src), destImage(fourier));
01881 
01882     // compute inverse Fourier transform
01883     // note that both source and destination image must be of type vigra::FFTWComplexImage
01884     vigra::FFTWComplexImage inverseFourier(w, h);
01885 
01886     fourierTransformInverse(srcImageRange(fourier), destImage(inverseFourier));
01887     \endcode
01888     
01889     new-style examples (using MultiArray):
01890     \code
01891     // compute Fourier transform of a real array, using the R2C algorithm
01892     MultiArray<2, double> src(Shape2(w, h));
01893     MultiArray<2, FFTWComplex<double> > fourier(fftwCorrespondingShapeR2C(src.shape()));
01894 
01895     fourierTransform(src, fourier);
01896 
01897     // compute inverse Fourier transform, using the C2R algorithm
01898     MultiArray<2, double> dest(src.shape());
01899     fourierTransformInverse(fourier, dest);
01900     \endcode
01901 
01902     \code
01903     // compute Fourier transform of a real array with standard algorithm
01904     MultiArray<2, double> src(Shape2(w, h));
01905     MultiArray<2, FFTWComplex<double> > fourier(src.shape());
01906 
01907     fourierTransform(src, fourier);
01908 
01909     // compute inverse Fourier transform, using the C2R algorithm
01910     MultiArray<2, double> dest(src.shape());
01911     fourierTransformInverse(fourier, dest);
01912     \endcode
01913     Complex input arrays are handled in the same way. 
01914     
01915 */
01916 doxygen_overloaded_function(template <...> void fourierTransform)
01917 
01918 inline void
01919 fourierTransform(FFTWComplexImage::const_traverser sul,
01920                  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01921                  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
01922 {
01923     detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
01924 }
01925 
01926 template <class SrcImageIterator, class SrcAccessor>
01927 void fourierTransform(SrcImageIterator srcUpperLeft,
01928                       SrcImageIterator srcLowerRight, SrcAccessor sa,
01929                       FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor da)
01930 {
01931     // copy real input images into a complex one...
01932     int w= srcLowerRight.x - srcUpperLeft.x;
01933     int h= srcLowerRight.y - srcUpperLeft.y;
01934 
01935     FFTWComplexImage workImage(w, h);
01936     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01937               destImage(workImage, FFTWWriteRealAccessor<>()));
01938 
01939     // ...and call the complex -> complex version of the algorithm
01940     FFTWComplexImage const & cworkImage = workImage;
01941     fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01942                      destUpperLeft, da);
01943 }
01944 
01945 template <class SrcImageIterator, class SrcAccessor>
01946 inline
01947 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01948                       pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
01949 {
01950     fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
01951 }
01952 
01953 /** \brief Compute inverse Fourier transforms.
01954 
01955     See \ref fourierTransform() for details.
01956 */
01957 doxygen_overloaded_function(template <...> void fourierTransformInverse)
01958 
01959 inline void
01960 fourierTransformInverse(FFTWComplexImage::const_traverser sul,
01961                         FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01962                         FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
01963 {
01964     detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
01965 }
01966 
01967 template <class DestImageIterator, class DestAccessor>
01968 void fourierTransformInverse(FFTWComplexImage::const_traverser sul,
01969                              FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01970                              DestImageIterator dul, DestAccessor dest)
01971 {
01972     int w = slr.x - sul.x;
01973     int h = slr.y - sul.y;
01974 
01975     FFTWComplexImage workImage(w, h);
01976     fourierTransformInverse(sul, slr, src, workImage.upperLeft(), workImage.accessor());
01977     copyImage(srcImageRange(workImage), destIter(dul, dest));
01978 }
01979 
01980 
01981 template <class DestImageIterator, class DestAccessor>
01982 inline void
01983 fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
01984                                FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
01985                         pair<DestImageIterator, DestAccessor> dest)
01986 {
01987     fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second);
01988 }
01989 
01990 
01991 
01992 /********************************************************/
01993 /*                                                      */
01994 /*                   applyFourierFilter                 */
01995 /*                                                      */
01996 /********************************************************/
01997 
01998 /** \brief Apply a filter (defined in the frequency domain) to an image.
01999 
02000     After transferring the image into the frequency domain, it is
02001     multiplied pixel-wise with the filter and transformed back. The
02002     result is put into the given destination image which must have the right size.
02003     The result will be normalized to compensate for the two FFTs.
02004 
02005     If the destination image is scalar, only the real part of the result image is
02006     retained. In this case, you are responsible for choosing a filter image
02007     which ensures a zero imaginary part of the result (e.g. use a real, even symmetric
02008     filter image, or a purely imaginary, odd symmetric one).
02009 
02010     The DC entry of the filter must be in the upper left, which is the
02011     position where FFTW expects it (see \ref moveDCToUpperLeft()).
02012     
02013     See also \ref convolveFFT() for corresponding functionality on the basis of the
02014     \ref MultiArrayView interface.
02015 
02016     <b> Declarations:</b>
02017 
02018     pass arguments explicitly:
02019     \code
02020     namespace vigra {
02021         template <class SrcImageIterator, class SrcAccessor,
02022                   class FilterImageIterator, class FilterAccessor,
02023                   class DestImageIterator, class DestAccessor>
02024         void applyFourierFilter(SrcImageIterator srcUpperLeft,
02025                                 SrcImageIterator srcLowerRight, SrcAccessor sa,
02026                                 FilterImageIterator filterUpperLeft, FilterAccessor fa,
02027                                 DestImageIterator destUpperLeft, DestAccessor da);
02028     }
02029     \endcode
02030 
02031     use argument objects in conjunction with \ref ArgumentObjectFactories :
02032     \code
02033     namespace vigra {
02034         template <class SrcImageIterator, class SrcAccessor,
02035                   class FilterImageIterator, class FilterAccessor,
02036                   class DestImageIterator, class DestAccessor>
02037         void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
02038                                 pair<FilterImageIterator, FilterAccessor> filter,
02039                                 pair<DestImageIterator, DestAccessor> dest);
02040     }
02041     \endcode
02042 
02043     <b> Usage:</b>
02044 
02045     <b>\#include</b> <vigra/fftw3.hxx><br>
02046     Namespace: vigra
02047 
02048     \code
02049     // create a Gaussian filter in Fourier space
02050     vigra::FImage gaussFilter(w, h), filter(w, h);
02051     for(int y=0; y<h; ++y)
02052         for(int x=0; x<w; ++x)
02053         {
02054             xx = float(x - w / 2) / w;
02055             yy = float(y - h / 2) / h;
02056 
02057             gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
02058         }
02059 
02060     // applyFourierFilter() expects the filter's DC in the upper left
02061     moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
02062 
02063     vigra::FFTWComplexImage result(w, h);
02064 
02065     vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
02066     \endcode
02067 
02068     For inspection of the result, \ref FFTWMagnitudeAccessor might be
02069     useful. If you want to apply the same filter repeatedly, it may be more
02070     efficient to use the FFTW functions directly with FFTW plans optimized
02071     for good performance.
02072 */
02073 doxygen_overloaded_function(template <...> void applyFourierFilter)
02074 
02075 template <class SrcImageIterator, class SrcAccessor,
02076           class FilterImageIterator, class FilterAccessor,
02077           class DestImageIterator, class DestAccessor>
02078 void applyFourierFilter(SrcImageIterator srcUpperLeft,
02079                         SrcImageIterator srcLowerRight, SrcAccessor sa,
02080                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
02081                         DestImageIterator destUpperLeft, DestAccessor da)
02082 {
02083     // copy real input images into a complex one...
02084     int w = int(srcLowerRight.x - srcUpperLeft.x);
02085     int h = int(srcLowerRight.y - srcUpperLeft.y);
02086 
02087     FFTWComplexImage workImage(w, h);
02088     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
02089               destImage(workImage, FFTWWriteRealAccessor<>()));
02090 
02091     // ...and call the impl
02092     FFTWComplexImage const & cworkImage = workImage;
02093     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
02094                            filterUpperLeft, fa,
02095                            destUpperLeft, da);
02096 }
02097 
02098 template <class FilterImageIterator, class FilterAccessor,
02099           class DestImageIterator, class DestAccessor>
02100 inline
02101 void applyFourierFilter(
02102     FFTWComplexImage::const_traverser srcUpperLeft,
02103     FFTWComplexImage::const_traverser srcLowerRight,
02104     FFTWComplexImage::ConstAccessor sa,
02105     FilterImageIterator filterUpperLeft, FilterAccessor fa,
02106     DestImageIterator destUpperLeft, DestAccessor da)
02107 {
02108     int w = srcLowerRight.x - srcUpperLeft.x;
02109     int h = srcLowerRight.y - srcUpperLeft.y;
02110 
02111     // test for right memory layout (fftw expects a 2*width*height floats array)
02112     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
02113         applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
02114                                filterUpperLeft, fa,
02115                                destUpperLeft, da);
02116     else
02117     {
02118         FFTWComplexImage workImage(w, h);
02119         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
02120                   destImage(workImage));
02121 
02122         FFTWComplexImage const & cworkImage = workImage;
02123         applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
02124                                filterUpperLeft, fa,
02125                                destUpperLeft, da);
02126     }
02127 }
02128 
02129 template <class SrcImageIterator, class SrcAccessor,
02130           class FilterImageIterator, class FilterAccessor,
02131           class DestImageIterator, class DestAccessor>
02132 inline
02133 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
02134                         pair<FilterImageIterator, FilterAccessor> filter,
02135                         pair<DestImageIterator, DestAccessor> dest)
02136 {
02137     applyFourierFilter(src.first, src.second, src.third,
02138                        filter.first, filter.second,
02139                        dest.first, dest.second);
02140 }
02141 
02142 template <class FilterImageIterator, class FilterAccessor,
02143           class DestImageIterator, class DestAccessor>
02144 void applyFourierFilterImpl(
02145     FFTWComplexImage::const_traverser srcUpperLeft,
02146     FFTWComplexImage::const_traverser srcLowerRight,
02147     FFTWComplexImage::ConstAccessor,
02148     FilterImageIterator filterUpperLeft, FilterAccessor fa,
02149     DestImageIterator destUpperLeft, DestAccessor da)
02150 {
02151     int w = int(srcLowerRight.x - srcUpperLeft.x);
02152     int h = int(srcLowerRight.y - srcUpperLeft.y);
02153 
02154     FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
02155 
02156     // FFT from srcImage to complexResultImg
02157     fftw_plan forwardPlan=
02158         fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
02159                                (fftw_complex *)complexResultImg.begin(),
02160                                FFTW_FORWARD, FFTW_ESTIMATE );
02161     fftw_execute(forwardPlan);
02162     fftw_destroy_plan(forwardPlan);
02163 
02164     // convolve in freq. domain (in complexResultImg)
02165     combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
02166                      destImage(complexResultImg), std::multiplies<FFTWComplex<> >());
02167 
02168     // FFT back into spatial domain (inplace in complexResultImg)
02169     fftw_plan backwardPlan=
02170         fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
02171                                (fftw_complex *)complexResultImg.begin(),
02172                                FFTW_BACKWARD, FFTW_ESTIMATE);
02173     fftw_execute(backwardPlan);
02174     fftw_destroy_plan(backwardPlan);
02175 
02176     typedef typename
02177         NumericTraits<typename DestAccessor::value_type>::isScalar
02178         isScalarResult;
02179 
02180     // normalization (after FFTs), maybe stripping imaginary part
02181     applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
02182                                         isScalarResult());
02183 }
02184 
02185 template <class DestImageIterator, class DestAccessor>
02186 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
02187                                          DestImageIterator destUpperLeft,
02188                                          DestAccessor da,
02189                                          VigraFalseType)
02190 {
02191     double normFactor= 1.0/(srcImage.width() * srcImage.height());
02192 
02193     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
02194     {
02195         DestImageIterator dIt= destUpperLeft;
02196         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
02197         {
02198             da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
02199             da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
02200         }
02201     }
02202 }
02203 
02204 inline
02205 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
02206         FFTWComplexImage::traverser destUpperLeft,
02207         FFTWComplexImage::Accessor da,
02208         VigraFalseType)
02209 {
02210     transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
02211                    linearIntensityTransform<FFTWComplex<> >(1.0/(srcImage.width() * srcImage.height())));
02212 }
02213 
02214 template <class DestImageIterator, class DestAccessor>
02215 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
02216                                          DestImageIterator destUpperLeft,
02217                                          DestAccessor da,
02218                                          VigraTrueType)
02219 {
02220     double normFactor= 1.0/(srcImage.width() * srcImage.height());
02221 
02222     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
02223     {
02224         DestImageIterator dIt= destUpperLeft;
02225         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
02226             da.set(srcImage(x, y).re()*normFactor, dIt);
02227     }
02228 }
02229 
02230 /**********************************************************/
02231 /*                                                        */
02232 /*                applyFourierFilterFamily                */
02233 /*                                                        */
02234 /**********************************************************/
02235 
02236 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
02237 
02238     This provides the same functionality as \ref applyFourierFilter(),
02239     but applying several filters at once allows to avoid
02240     repeated Fourier transforms of the source image.
02241 
02242     Filters and result images must be stored in \ref vigra::ImageArray data
02243     structures. In contrast to \ref applyFourierFilter(), this function adjusts
02244     the size of the result images and the the length of the array.
02245 
02246     <b> Declarations:</b>
02247 
02248     pass arguments explicitly:
02249     \code
02250     namespace vigra {
02251         template <class SrcImageIterator, class SrcAccessor, class FilterType>
02252         void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
02253                                       SrcImageIterator srcLowerRight, SrcAccessor sa,
02254                                       const ImageArray<FilterType> &filters,
02255                                       ImageArray<FFTWComplexImage> &results)
02256     }
02257     \endcode
02258 
02259     use argument objects in conjunction with \ref ArgumentObjectFactories :
02260     \code
02261     namespace vigra {
02262         template <class SrcImageIterator, class SrcAccessor, class FilterType>
02263         void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
02264                                       const ImageArray<FilterType> &filters,
02265                                       ImageArray<FFTWComplexImage> &results)
02266     }
02267     \endcode
02268 
02269     <b> Usage:</b>
02270 
02271     <b>\#include</b> <vigra/fftw3.hxx><br>
02272     Namespace: vigra
02273 
02274     \code
02275     // assuming the presence of a real-valued image named "image" to
02276     // be filtered in this example
02277 
02278     vigra::ImageArray<vigra::FImage> filters(16, image.size());
02279 
02280     for (int i=0; i<filters.size(); i++)
02281          // create some meaningful filters here
02282          createMyFilterOfScale(i, destImage(filters[i]));
02283 
02284     vigra::ImageArray<vigra::FFTWComplexImage> results();
02285 
02286     vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
02287     \endcode
02288 */
02289 doxygen_overloaded_function(template <...> void applyFourierFilterFamily)
02290 
02291 template <class SrcImageIterator, class SrcAccessor,
02292           class FilterType, class DestImage>
02293 inline
02294 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
02295                               const ImageArray<FilterType> &filters,
02296                               ImageArray<DestImage> &results)
02297 {
02298     applyFourierFilterFamily(src.first, src.second, src.third,
02299                              filters, results);
02300 }
02301 
02302 template <class SrcImageIterator, class SrcAccessor,
02303           class FilterType, class DestImage>
02304 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
02305                               SrcImageIterator srcLowerRight, SrcAccessor sa,
02306                               const ImageArray<FilterType> &filters,
02307                               ImageArray<DestImage> &results)
02308 {
02309     int w = int(srcLowerRight.x - srcUpperLeft.x);
02310     int h = int(srcLowerRight.y - srcUpperLeft.y);
02311 
02312     FFTWComplexImage workImage(w, h);
02313     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
02314               destImage(workImage, FFTWWriteRealAccessor<>()));
02315 
02316     FFTWComplexImage const & cworkImage = workImage;
02317     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
02318                                  filters, results);
02319 }
02320 
02321 template <class FilterType, class DestImage>
02322 inline
02323 void applyFourierFilterFamily(
02324     FFTWComplexImage::const_traverser srcUpperLeft,
02325     FFTWComplexImage::const_traverser srcLowerRight,
02326     FFTWComplexImage::ConstAccessor sa,
02327     const ImageArray<FilterType> &filters,
02328     ImageArray<DestImage> &results)
02329 {
02330     int w= srcLowerRight.x - srcUpperLeft.x;
02331 
02332     // test for right memory layout (fftw expects a 2*width*height floats array)
02333     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
02334         applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
02335                                      filters, results);
02336     else
02337     {
02338         int h = srcLowerRight.y - srcUpperLeft.y;
02339         FFTWComplexImage workImage(w, h);
02340         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
02341                   destImage(workImage));
02342 
02343         FFTWComplexImage const & cworkImage = workImage;
02344         applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
02345                                      filters, results);
02346     }
02347 }
02348 
02349 template <class FilterType, class DestImage>
02350 void applyFourierFilterFamilyImpl(
02351     FFTWComplexImage::const_traverser srcUpperLeft,
02352     FFTWComplexImage::const_traverser srcLowerRight,
02353     FFTWComplexImage::ConstAccessor sa,
02354     const ImageArray<FilterType> &filters,
02355     ImageArray<DestImage> &results)
02356 {
02357     // FIXME: sa is not used
02358     // (maybe check if StandardAccessor, else copy?)    
02359 
02360     // make sure the filter images have the right dimensions
02361     vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
02362                        "applyFourierFilterFamily called with src image size != filters.imageSize()!");
02363 
02364     // make sure the result image array has the right dimensions
02365     results.resize(filters.size());
02366     results.resizeImages(filters.imageSize());
02367 
02368     // FFT from srcImage to freqImage
02369     int w = int(srcLowerRight.x - srcUpperLeft.x);
02370     int h = int(srcLowerRight.y - srcUpperLeft.y);
02371 
02372     FFTWComplexImage freqImage(w, h);
02373     FFTWComplexImage result(w, h);
02374 
02375     fftw_plan forwardPlan=
02376         fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
02377                                (fftw_complex *)freqImage.begin(),
02378                                FFTW_FORWARD, FFTW_ESTIMATE );
02379     fftw_execute(forwardPlan);
02380     fftw_destroy_plan(forwardPlan);
02381 
02382     fftw_plan backwardPlan=
02383         fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
02384                                (fftw_complex *)result.begin(),
02385                                FFTW_BACKWARD, FFTW_ESTIMATE );
02386     typedef typename
02387         NumericTraits<typename DestImage::Accessor::value_type>::isScalar
02388         isScalarResult;
02389 
02390     // convolve with filters in freq. domain
02391     for (unsigned int i= 0;  i < filters.size(); i++)
02392     {
02393         combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
02394                          destImage(result), std::multiplies<FFTWComplex<> >());
02395 
02396         // FFT back into spatial domain (inplace in destImage)
02397         fftw_execute(backwardPlan);
02398 
02399         // normalization (after FFTs), maybe stripping imaginary part
02400         applyFourierFilterImplNormalization(result,
02401                                             results[i].upperLeft(), results[i].accessor(),
02402                                             isScalarResult());
02403     }
02404     fftw_destroy_plan(backwardPlan);
02405 }
02406 
02407 /********************************************************/
02408 /*                                                      */
02409 /*                fourierTransformReal                  */
02410 /*                                                      */
02411 /********************************************************/
02412 
02413 /** \brief Real Fourier transforms for even and odd boundary conditions
02414            (aka. cosine and sine transforms).
02415 
02416 
02417     If the image is real and has even symmetry, its Fourier transform
02418     is also real and has even symmetry. The Fourier transform of a real image with odd
02419     symmetry is imaginary and has odd symmetry. In either case, only about a quarter
02420     of the pixels need to be stored because the rest can be calculated from the symmetry
02421     properties. This is especially useful, if the original image is implicitly assumed
02422     to have reflective or anti-reflective boundary conditions. Then the "negative"
02423     pixel locations are defined as
02424 
02425     \code
02426     even (reflective boundary conditions):      f[-x] = f[x]     (x = 1,...,N-1)
02427     odd (anti-reflective boundary conditions):  f[-1] = 0
02428                                                 f[-x] = -f[x-2]  (x = 2,...,N-1)
02429     \endcode
02430 
02431     end similar at the other boundary (see the FFTW documentation for details).
02432     This has the advantage that more efficient Fourier transforms that use only
02433     real numbers can be implemented. These are also known as cosine and sine transforms
02434     respectively.
02435 
02436     If you use the odd transform it is important to note that in the Fourier domain,
02437     the DC component is always zero and is therefore dropped from the data structure.
02438     This means that index 0 in an odd symmetric Fourier domain image refers to
02439     the <i>first</i> harmonic. This is especially important if an image is first
02440     cosine transformed (even symmetry), then in the Fourier domain multiplied
02441     with an odd symmetric filter (e.g. a first derivative) and finally transformed
02442     back to the spatial domain with a sine transform (odd symmetric). For this to work
02443     properly the image must be shifted left or up by one pixel (depending on whether
02444     the x- or y-axis is odd symmetric) before the inverse transform can be applied.
02445     (see example below).
02446 
02447     The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
02448     where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
02449     whether the x- and y-axis is to be transformed using even or odd symmetry.
02450     The same functions can be used for both the forward and inverse transforms,
02451     only the normalization changes. For signal processing, the following
02452     normalization factors are most appropriate:
02453 
02454     \code
02455                           forward             inverse
02456     ------------------------------------------------------------
02457     X even, Y even           1.0         4.0 * (w-1) * (h-1)
02458     X even, Y odd           -1.0        -4.0 * (w-1) * (h+1)
02459     X odd,  Y even          -1.0        -4.0 * (w+1) * (h-1)
02460     X odd,  Y odd            1.0         4.0 * (w+1) * (h+1)
02461     \endcode
02462 
02463     where <tt>w</tt> and <tt>h</tt> denote the image width and height.
02464 
02465     <b> Declarations:</b>
02466 
02467     pass arguments explicitly:
02468     \code
02469     namespace vigra {
02470         template <class SrcTraverser, class SrcAccessor,
02471                   class DestTraverser, class DestAccessor>
02472         void
02473         fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
02474                                DestTraverser dul, DestAccessor dest, fftw_real norm);
02475 
02476         fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
02477     }
02478     \endcode
02479 
02480 
02481     use argument objects in conjunction with \ref ArgumentObjectFactories :
02482     \code
02483     namespace vigra {
02484         template <class SrcTraverser, class SrcAccessor,
02485                   class DestTraverser, class DestAccessor>
02486         void
02487         fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
02488                                pair<DestTraverser, DestAccessor> dest, fftw_real norm);
02489 
02490         fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
02491     }
02492     \endcode
02493 
02494     <b> Usage:</b>
02495 
02496         <b>\#include</b> <vigra/fftw3.hxx><br>
02497         Namespace: vigra
02498 
02499     \code
02500     vigra::FImage spatial(width,height), fourier(width,height);
02501     ... // fill image with data
02502 
02503     // forward cosine transform == reflective boundary conditions
02504     fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
02505 
02506     // multiply with a first derivative of Gaussian in x-direction
02507     for(int y = 0; y < height; ++y)
02508     {
02509         for(int x = 1; x < width; ++x)
02510         {
02511             double dx = x * M_PI / (width - 1);
02512             double dy = y * M_PI / (height - 1);
02513             fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
02514         }
02515         fourier(width-1, y) = 0.0;
02516     }
02517 
02518     // inverse transform -- odd symmetry in x-direction, even in y,
02519     //                      due to symmetry of the filter
02520     fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
02521                            (fftw_real)-4.0 * (width+1) * (height-1));
02522     \endcode
02523 */
02524 doxygen_overloaded_function(template <...> void fourierTransformReal)
02525 
02526 template <class SrcTraverser, class SrcAccessor,
02527           class DestTraverser, class DestAccessor>
02528 inline void
02529 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
02530                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
02531 {
02532     fourierTransformRealEE(src.first, src.second, src.third,
02533                                    dest.first, dest.second, norm);
02534 }
02535 
02536 template <class SrcTraverser, class SrcAccessor,
02537           class DestTraverser, class DestAccessor>
02538 inline void
02539 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
02540                                DestTraverser dul, DestAccessor dest, fftw_real norm)
02541 {
02542     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02543                                       norm, FFTW_REDFT00, FFTW_REDFT00);
02544 }
02545 
02546 template <class DestTraverser, class DestAccessor>
02547 inline void
02548 fourierTransformRealEE(
02549          FFTWRealImage::const_traverser sul,
02550          FFTWRealImage::const_traverser slr,
02551          FFTWRealImage::Accessor src,
02552          DestTraverser dul, DestAccessor dest, fftw_real norm)
02553 {
02554     int w = slr.x - sul.x;
02555 
02556     // test for right memory layout (fftw expects a width*height fftw_real array)
02557     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
02558         fourierTransformRealImpl(sul, slr, dul, dest,
02559                                  norm, FFTW_REDFT00, FFTW_REDFT00);
02560     else
02561         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02562                                  norm, FFTW_REDFT00, FFTW_REDFT00);
02563 }
02564 
02565 /********************************************************************/
02566 
02567 template <class SrcTraverser, class SrcAccessor,
02568           class DestTraverser, class DestAccessor>
02569 inline void
02570 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
02571                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
02572 {
02573     fourierTransformRealOE(src.first, src.second, src.third,
02574                                    dest.first, dest.second, norm);
02575 }
02576 
02577 template <class SrcTraverser, class SrcAccessor,
02578           class DestTraverser, class DestAccessor>
02579 inline void
02580 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
02581                                DestTraverser dul, DestAccessor dest, fftw_real norm)
02582 {
02583     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02584                                       norm, FFTW_RODFT00, FFTW_REDFT00);
02585 }
02586 
02587 template <class DestTraverser, class DestAccessor>
02588 inline void
02589 fourierTransformRealOE(
02590          FFTWRealImage::const_traverser sul,
02591          FFTWRealImage::const_traverser slr,
02592          FFTWRealImage::Accessor src,
02593          DestTraverser dul, DestAccessor dest, fftw_real norm)
02594 {
02595     int w = slr.x - sul.x;
02596 
02597     // test for right memory layout (fftw expects a width*height fftw_real array)
02598     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
02599         fourierTransformRealImpl(sul, slr, dul, dest,
02600                                  norm, FFTW_RODFT00, FFTW_REDFT00);
02601     else
02602         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02603                                  norm, FFTW_RODFT00, FFTW_REDFT00);
02604 }
02605 
02606 /********************************************************************/
02607 
02608 template <class SrcTraverser, class SrcAccessor,
02609           class DestTraverser, class DestAccessor>
02610 inline void
02611 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
02612                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
02613 {
02614     fourierTransformRealEO(src.first, src.second, src.third,
02615                                    dest.first, dest.second, norm);
02616 }
02617 
02618 template <class SrcTraverser, class SrcAccessor,
02619           class DestTraverser, class DestAccessor>
02620 inline void
02621 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
02622                                DestTraverser dul, DestAccessor dest, fftw_real norm)
02623 {
02624     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02625                                       norm, FFTW_REDFT00, FFTW_RODFT00);
02626 }
02627 
02628 template <class DestTraverser, class DestAccessor>
02629 inline void
02630 fourierTransformRealEO(
02631          FFTWRealImage::const_traverser sul,
02632          FFTWRealImage::const_traverser slr,
02633          FFTWRealImage::Accessor src,
02634          DestTraverser dul, DestAccessor dest, fftw_real norm)
02635 {
02636     int w = slr.x - sul.x;
02637 
02638     // test for right memory layout (fftw expects a width*height fftw_real array)
02639     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
02640         fourierTransformRealImpl(sul, slr, dul, dest,
02641                                  norm, FFTW_REDFT00, FFTW_RODFT00);
02642     else
02643         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02644                                  norm, FFTW_REDFT00, FFTW_RODFT00);
02645 }
02646 
02647 /********************************************************************/
02648 
02649 template <class SrcTraverser, class SrcAccessor,
02650           class DestTraverser, class DestAccessor>
02651 inline void
02652 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
02653                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
02654 {
02655     fourierTransformRealOO(src.first, src.second, src.third,
02656                                    dest.first, dest.second, norm);
02657 }
02658 
02659 template <class SrcTraverser, class SrcAccessor,
02660           class DestTraverser, class DestAccessor>
02661 inline void
02662 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
02663                                DestTraverser dul, DestAccessor dest, fftw_real norm)
02664 {
02665     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02666                                       norm, FFTW_RODFT00, FFTW_RODFT00);
02667 }
02668 
02669 template <class DestTraverser, class DestAccessor>
02670 inline void
02671 fourierTransformRealOO(
02672          FFTWRealImage::const_traverser sul,
02673          FFTWRealImage::const_traverser slr,
02674          FFTWRealImage::Accessor src,
02675          DestTraverser dul, DestAccessor dest, fftw_real norm)
02676 {
02677     int w = slr.x - sul.x;
02678 
02679     // test for right memory layout (fftw expects a width*height fftw_real array)
02680     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
02681         fourierTransformRealImpl(sul, slr, dul, dest,
02682                                  norm, FFTW_RODFT00, FFTW_RODFT00);
02683     else
02684         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
02685                                  norm, FFTW_RODFT00, FFTW_RODFT00);
02686 }
02687 
02688 /*******************************************************************/
02689 
02690 template <class SrcTraverser, class SrcAccessor,
02691           class DestTraverser, class DestAccessor>
02692 void
02693 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
02694                                   DestTraverser dul, DestAccessor dest,
02695                                   fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
02696 {
02697     FFTWRealImage workImage(slr - sul);
02698     copyImage(srcIterRange(sul, slr, src), destImage(workImage));
02699     FFTWRealImage const & cworkImage = workImage;
02700     fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
02701                              dul, dest, norm, kindx, kindy);
02702 }
02703 
02704 
02705 template <class DestTraverser, class DestAccessor>
02706 void
02707 fourierTransformRealImpl(
02708          FFTWRealImage::const_traverser sul,
02709          FFTWRealImage::const_traverser slr,
02710          DestTraverser dul, DestAccessor dest,
02711          fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
02712 {
02713     int w = slr.x - sul.x;
02714     int h = slr.y - sul.y;
02715     BasicImage<fftw_real> res(w, h);
02716 
02717     fftw_plan plan = fftw_plan_r2r_2d(h, w,
02718                          (fftw_real *)&(*sul), (fftw_real *)res.begin(),
02719                          kindy, kindx, FFTW_ESTIMATE);
02720     fftw_execute(plan);
02721     fftw_destroy_plan(plan);
02722 
02723     if(norm != 1.0)
02724         transformImage(srcImageRange(res), destIter(dul, dest),
02725                        std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
02726     else
02727         copyImage(srcImageRange(res), destIter(dul, dest));
02728 }
02729 
02730 
02731 //@}
02732 
02733 } // namespace vigra
02734 
02735 #endif // VIGRA_FFTW3_HXX

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

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)