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

vigra/fftw3.hxx |

    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)