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

| vigra/fftw3.hxx |
--- | --- | ---

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.3.2, Jan 27 2005 )                                    */
    00008 /*    You may use, modify, and distribute this software according       */
    00009 /*    to the terms stated in the LICENSE file included in               */
    00010 /*    the VIGRA distribution.                                           */
    00011 /*                                                                      */
    00012 /*    The VIGRA Website is                                              */
    00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00014 /*    Please direct questions, bug reports, and contributions to        */
    00015 /*        koethe@informatik.uni-hamburg.de                              */
    00016 /*                                                                      */
    00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
    00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
    00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
    00020 /*                                                                      */
    00021 /************************************************************************/
    00022
    00023 #ifndef VIGRA_FFTW3_HXX
    00024 #define VIGRA_FFTW3_HXX
    00025
    00026 #include <cmath>
    00027 #include <functional>
    00028 #include "vigra/stdimage.hxx"
    00029 #include "vigra/copyimage.hxx"
    00030 #include "vigra/transformimage.hxx"
    00031 #include "vigra/combineimages.hxx"
    00032 #include "vigra/numerictraits.hxx"
    00033 #include "vigra/imagecontainer.hxx"
    00034 #include <fftw3.h>
    00035
    00036 namespace vigra {
    00037
    00038 typedef double fftw_real;
    00039
    00040 /********************************************************/
    00041 /*                                                      */
    00042 /*                    FFTWComplex                       */
    00043 /*                                                      */
    00044 /********************************************************/
    00045
    00046 /** \brief Wrapper class for the FFTW type '<TT>fftw_complex</TT>'.
    00047
    00048     This class provides constructors and other member functions
    00049     for the C struct '<TT>fftw_complex</TT>'. This struct is the basic
    00050     pixel type of the <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a>
    00051     library. It inherits the data members '<TT>re</TT>' and '<TT>im</TT>'
    00052     that denote the real and imaginary part of the number. In addition it
    00053     defines transformations to polar coordinates,
    00054     as well as \ref FFTWComplexOperators "arithmetic operators" and
    00055     \ref FFTWComplexAccessors "accessors".
    00056
    00057     FFTWComplex implements the concepts \ref AlgebraicField and
    00058     \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt>
    00059     and <tt>FFTWComplexImage</tt> are defined.
    00060
    00061     <b>See also:</b>
    00062     <ul>
    00063         <li> \ref FFTWComplexTraits
    00064         <li> \ref FFTWComplexOperators
    00065         <li> \ref FFTWComplexAccessors
    00066     </ul>
    00067
    00068     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00069     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00070     Namespace: vigra
    00071 */
    00072 class FFTWComplex
    00073 {
    00074     fftw_complex data_;
    00075
    00076   public:
    00077         /** The complex' component type, as defined in '<TT>fftw3.h</TT>'
    00078         */
    00079     typedef fftw_real value_type;
    00080
    00081         /** reference type (result of operator[])
    00082         */
    00083     typedef fftw_real & reference;
    00084
    00085         /** const reference type (result of operator[] const)
    00086         */
    00087     typedef fftw_real const & const_reference;
    00088
    00089         /** const reference type (result of operator[] const)
    00090         */
    00091     typedef fftw_real * iterator;
    00092
    00093         /** const reference type (result of operator[] const)
    00094         */
    00095     typedef fftw_real const * const_iterator;
    00096
    00097         /** Construct from real and imaginary part.
    00098             Default: 0.
    00099         */
    00100     FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0)
    00101     {
    00102         data_[0] = re;
    00103         data_[1] = im;
    00104     }
    00105
    00106         /** Copy constructor.
    00107         */
    00108     FFTWComplex(FFTWComplex const & o)
    00109     {
    00110         data_[0] = o.data_[0];
    00111         data_[1] = o.data_[1];
    00112     }
    00113
    00114         /** Construct from plain <TT>fftw_complex</TT>.
    00115         */
    00116     FFTWComplex(fftw_complex const & o)
    00117     {
    00118         data_[0] = o[0];
    00119         data_[1] = o[1];
    00120     }
    00121
    00122         /** Construct from TinyVector.
    00123         */
    00124     template <class T>
    00125     FFTWComplex(TinyVector<T, 2> const & o)
    00126     {
    00127         data_[0] = o[0];
    00128         data_[1] = o[1];
    00129     }
    00130
    00131         /** Assignment.
    00132         */
    00133     FFTWComplex& operator=(FFTWComplex const & o)
    00134     {
    00135         data_[0] = o.data_[0];
    00136         data_[1] = o.data_[1];
    00137         return *this;
    00138     }
    00139
    00140         /** Assignment.
    00141         */
    00142     FFTWComplex& operator=(fftw_complex const & o)
    00143     {
    00144         data_[0] = o[0];
    00145         data_[1] = o[1];
    00146         return *this;
    00147     }
    00148
    00149         /** Assignment.
    00150         */
    00151     FFTWComplex& operator=(fftw_real const & o)
    00152     {
    00153         data_[0] = o;
    00154         data_[1] = 0.0;
    00155         return *this;
    00156     }
    00157
    00158         /** Assignment.
    00159         */
    00160     template <class T>
    00161     FFTWComplex& operator=(TinyVector<T, 2> const & o)
    00162     {
    00163         data_[0] = o[0];
    00164         data_[1] = o[1];
    00165         return *this;
    00166     }
    00167
    00168     reference re()
    00169         { return data_[0]; }
    00170
    00171     const_reference re() const
    00172         { return data_[0]; }
    00173
    00174     reference im()
    00175         { return data_[1]; }
    00176
    00177     const_reference im() const
    00178         { return data_[1]; }
    00179
    00180         /** Unary negation.
    00181         */
    00182     FFTWComplex operator-() const
    00183         { return FFTWComplex(-data_[0], -data_[1]); }
    00184
    00185         /** Squared magnitude x*conj(x)
    00186         */
    00187     value_type squaredMagnitude() const
    00188         { return data_[0]*data_[0]+data_[1]*data_[1]; }
    00189
    00190         /** Magnitude (length of radius vector).
    00191         */
    00192     value_type magnitude() const
    00193         { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
    00194
    00195         /** Phase angle.
    00196         */
    00197     value_type phase() const
    00198         { return VIGRA_CSTD::atan2(data_[1], data_[0]); }
    00199
    00200         /** Access components as if number were a vector.
    00201         */
    00202     reference operator[](int i)
    00203         { return data_[i]; }
    00204
    00205         /** Read components as if number were a vector.
    00206         */
    00207     const_reference operator[](int i) const
    00208         { return data_[i]; }
    00209
    00210         /** Length of complex number (always 2).
    00211         */
    00212     int size() const
    00213         { return 2; }
    00214
    00215     iterator begin()
    00216         { return data_; }
    00217
    00218     iterator end()
    00219         { return data_ + 2; }
    00220
    00221     const_iterator begin() const
    00222         { return data_; }
    00223
    00224     const_iterator end() const
    00225         { return data_ + 2; }
    00226 };
    00227
    00228 /********************************************************/
    00229 /*                                                      */
    00230 /*                     FFTWComplexTraits                */
    00231 /*                                                      */
    00232 /********************************************************/
    00233
    00234 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex
    00235
    00236     The numeric and promote traits for fftw_complex and FFTWComplex follow
    00237     the general specifications for \ref NumericPromotionTraits and
    00238     \ref AlgebraicField. They are explicitly specialized for the types
    00239     involved:
    00240
    00241     \code
    00242
    00243     template<>
    00244     struct NumericTraits<fftw_complex>
    00245     {
    00246         typedef fftw_complex Promote;
    00247         typedef fftw_complex RealPromote;
    00248         typedef fftw_complex ComplexPromote;
    00249         typedef fftw_real    ValueType;
    00250
    00251         typedef VigraFalseType isIntegral;
    00252         typedef VigraFalseType isScalar;
    00253         typedef VigraFalseType isOrdered;
    00254         typedef VigraTrueType  isComplex;
    00255
    00256         // etc.
    00257     };
    00258
    00259     template<>
    00260     struct NumericTraits<FFTWComplex>
    00261     {
    00262         typedef FFTWComplex Promote;
    00263         typedef FFTWComplex RealPromote;
    00264         typedef FFTWComplex ComplexPromote;
    00265         typedef fftw_real   ValueType;
    00266
    00267         typedef VigraFalseType isIntegral;
    00268         typedef VigraFalseType isScalar;
    00269         typedef VigraFalseType isOrdered;
    00270         typedef VigraTrueType  isComplex;
    00271
    00272         // etc.
    00273     };
    00274
    00275     template <>
    00276     struct PromoteTraits<fftw_complex, fftw_complex>
    00277     {
    00278         typedef fftw_complex Promote;
    00279     };
    00280
    00281     template <>
    00282     struct PromoteTraits<fftw_complex, double>
    00283     {
    00284         typedef fftw_complex Promote;
    00285     };
    00286
    00287     template <>
    00288     struct PromoteTraits<double, fftw_complex>
    00289     {
    00290         typedef fftw_complex Promote;
    00291     };
    00292
    00293     template <>
    00294     struct PromoteTraits<FFTWComplex, FFTWComplex>
    00295     {
    00296         typedef FFTWComplex Promote;
    00297     };
    00298
    00299     template <>
    00300     struct PromoteTraits<FFTWComplex, double>
    00301     {
    00302         typedef FFTWComplex Promote;
    00303     };
    00304
    00305     template <>
    00306     struct PromoteTraits<double, FFTWComplex>
    00307     {
    00308         typedef FFTWComplex Promote;
    00309     };
    00310     \endcode
    00311
    00312     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00313     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00314     Namespace: vigra
    00315
    00316 */
    00317 template<>
    00318 struct NumericTraits<fftw_complex>
    00319 {
    00320     typedef fftw_complex Type;
    00321     typedef fftw_complex Promote;
    00322     typedef fftw_complex RealPromote;
    00323     typedef fftw_complex ComplexPromote;
    00324     typedef fftw_real    ValueType;
    00325     typedef VigraFalseType isIntegral;
    00326     typedef VigraFalseType isScalar;
    00327     typedef VigraFalseType isOrdered;
    00328     typedef VigraTrueType  isComplex;
    00329
    00330     static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
    00331     static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
    00332     static FFTWComplex nonZero() { return one(); }
    00333
    00334     static const Promote & toPromote(const Type & v) { return v; }
    00335     static const RealPromote & toRealPromote(const Type & v) { return v; }
    00336     static const Type & fromPromote(const Promote & v) { return v; }
    00337     static const Type & fromRealPromote(const RealPromote & v) { return v; }
    00338 };
    00339
    00340 template<>
    00341 struct NumericTraits<FFTWComplex>
    00342 {
    00343     typedef FFTWComplex Type;
    00344     typedef FFTWComplex Promote;
    00345     typedef FFTWComplex RealPromote;
    00346     typedef FFTWComplex ComplexPromote;
    00347     typedef fftw_real   ValueType;
    00348     typedef VigraFalseType isIntegral;
    00349     typedef VigraFalseType isScalar;
    00350     typedef VigraFalseType isOrdered;
    00351     typedef VigraTrueType  isComplex;
    00352
    00353     static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
    00354     static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
    00355     static FFTWComplex nonZero() { return one(); }
    00356
    00357     static const Promote & toPromote(const Type & v) { return v; }
    00358     static const RealPromote & toRealPromote(const Type & v) { return v; }
    00359     static const Type & fromPromote(const Promote & v) { return v; }
    00360     static const Type & fromRealPromote(const RealPromote & v) { return v; }
    00361 };
    00362
    00363 template <>
    00364 struct PromoteTraits<fftw_complex, fftw_complex>
    00365 {
    00366     typedef fftw_complex Promote;
    00367 };
    00368
    00369 template <>
    00370 struct PromoteTraits<fftw_complex, double>
    00371 {
    00372     typedef fftw_complex Promote;
    00373 };
    00374
    00375 template <>
    00376 struct PromoteTraits<double, fftw_complex>
    00377 {
    00378     typedef fftw_complex Promote;
    00379 };
    00380
    00381 template <>
    00382 struct PromoteTraits<FFTWComplex, FFTWComplex>
    00383 {
    00384     typedef FFTWComplex Promote;
    00385 };
    00386
    00387 template <>
    00388 struct PromoteTraits<FFTWComplex, double>
    00389 {
    00390     typedef FFTWComplex Promote;
    00391 };
    00392
    00393 template <>
    00394 struct PromoteTraits<double, FFTWComplex>
    00395 {
    00396     typedef FFTWComplex Promote;
    00397 };
    00398
    00399
    00400 /********************************************************/
    00401 /*                                                      */
    00402 /*                    FFTWComplex Operations            */
    00403 /*                                                      */
    00404 /********************************************************/
    00405
    00406 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex
    00407
    00408     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00409     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00410
    00411     These functions fulfill the requirements of an Algebraic Field.
    00412     Return types are determined according to \ref FFTWComplexTraits.
    00413
    00414     Namespace: vigra
    00415     <p>
    00416
    00417  */
    00418 //@{
    00419     /// equal
    00420 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) {
    00421     return a.re() == b.re() && a.im() == b.im();
    00422 }
    00423
    00424     /// not equal
    00425 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) {
    00426     return a.re() != b.re() || a.im() != b.im();
    00427 }
    00428
    00429     /// add-assignment
    00430 inline FFTWComplex &operator +=(FFTWComplex &a, const FFTWComplex &b) {
    00431     a.re() += b.re();
    00432     a.im() += b.im();
    00433     return a;
    00434 }
    00435
    00436     /// subtract-assignment
    00437 inline FFTWComplex &operator -=(FFTWComplex &a, const FFTWComplex &b) {
    00438     a.re() -= b.re();
    00439     a.im() -= b.im();
    00440     return a;
    00441 }
    00442
    00443     /// multiply-assignment
    00444 inline FFTWComplex &operator *=(FFTWComplex &a, const FFTWComplex &b) {
    00445     FFTWComplex::value_type t = a.re()*b.re()-a.im()*b.im();
    00446     a.im() = a.re()*b.im()+a.im()*b.re();
    00447     a.re() = t;
    00448     return a;
    00449 }
    00450
    00451     /// divide-assignment
    00452 inline FFTWComplex &operator /=(FFTWComplex &a, const FFTWComplex &b) {
    00453     FFTWComplex::value_type sm = b.squaredMagnitude();
    00454     FFTWComplex::value_type t = (a.re()*b.re()+a.im()*b.im())/sm;
    00455     a.im() = (b.re()*a.im()-a.re()*b.im())/sm;
    00456     a.re() = t;
    00457     return a;
    00458 }
    00459
    00460     /// multiply-assignment with scalar double
    00461 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) {
    00462     a.re() *= b;
    00463     a.im() *= b;
    00464     return a;
    00465 }
    00466
    00467     /// divide-assignment with scalar double
    00468 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) {
    00469     a.re() /= b;
    00470     a.im() /= b;
    00471     return a;
    00472 }
    00473
    00474     /// addition
    00475 inline FFTWComplex operator +(FFTWComplex a, const FFTWComplex &b) {
    00476     a += b;
    00477     return a;
    00478 }
    00479
    00480     /// subtraction
    00481 inline FFTWComplex operator -(FFTWComplex a, const FFTWComplex &b) {
    00482     a -= b;
    00483     return a;
    00484 }
    00485
    00486     /// multiplication
    00487 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) {
    00488     a *= b;
    00489     return a;
    00490 }
    00491
    00492     /// right multiplication with scalar double
    00493 inline FFTWComplex operator *(FFTWComplex a, const double &b) {
    00494     a *= b;
    00495     return a;
    00496 }
    00497
    00498     /// left multiplication with scalar double
    00499 inline FFTWComplex operator *(const double &a, FFTWComplex b) {
    00500     b *= a;
    00501     return b;
    00502 }
    00503
    00504     /// division
    00505 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) {
    00506     a /= b;
    00507     return a;
    00508 }
    00509
    00510     /// right division with scalar double
    00511 inline FFTWComplex operator /(FFTWComplex a, const double &b) {
    00512     a /= b;
    00513     return a;
    00514 }
    00515
    00516 using VIGRA_CSTD::abs;
    00517
    00518     /// absolute value (= magnitude)
    00519 inline FFTWComplex::value_type abs(const FFTWComplex &a)
    00520 {
    00521     return a.magnitude();
    00522 }
    00523
    00524     /// complex conjugate
    00525 inline FFTWComplex conj(const FFTWComplex &a)
    00526 {
    00527     return FFTWComplex(a.re(), -a.im());
    00528 }
    00529
    00530 //@}
    00531
    00532
    00533 /** \addtogroup StandardImageTypes
    00534 */
    00535 //@{
    00536
    00537 /********************************************************/
    00538 /*                                                      */
    00539 /*                      FFTWRealImage                   */
    00540 /*                                                      */
    00541 /********************************************************/
    00542
    00543     /** Float (<tt>fftw_real</tt>) image.
    00544
    00545         The type <tt>fftw_real</tt> is defined as <tt>double</tt> (in FFTW 2 it used to be
    00546         either <tt>float</tt> or <tt>double</tt>, as specified during compilation of FFTW).
    00547         FFTWRealImage uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
    00548         their const counterparts to access the data.
    00549
    00550         <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00551         <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00552         Namespace: vigra
    00553     */
    00554 typedef BasicImage<fftw_real> FFTWRealImage;
    00555
    00556 /********************************************************/
    00557 /*                                                      */
    00558 /*                     FFTWComplexImage                 */
    00559 /*                                                      */
    00560 /********************************************************/
    00561
    00562 template<>
    00563 struct IteratorTraits<
    00564         BasicImageIterator<FFTWComplex, FFTWComplex **> >
    00565 {
    00566     typedef BasicImageIterator<FFTWComplex, FFTWComplex **>  Iterator;
    00567     typedef Iterator                             iterator;
    00568     typedef iterator::iterator_category          iterator_category;
    00569     typedef iterator::value_type                 value_type;
    00570     typedef iterator::reference                  reference;
    00571     typedef iterator::index_reference            index_reference;
    00572     typedef iterator::pointer                    pointer;
    00573     typedef iterator::difference_type            difference_type;
    00574     typedef iterator::row_iterator               row_iterator;
    00575     typedef iterator::column_iterator            column_iterator;
    00576     typedef VectorAccessor<FFTWComplex>          default_accessor;
    00577     typedef VectorAccessor<FFTWComplex>          DefaultAccessor;
    00578 };
    00579
    00580 template<>
    00581 struct IteratorTraits<
    00582         ConstBasicImageIterator<FFTWComplex, FFTWComplex **> >
    00583 {
    00584     typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **>    Iterator;
    00585     typedef Iterator                             iterator;
    00586     typedef iterator::iterator_category          iterator_category;
    00587     typedef iterator::value_type                 value_type;
    00588     typedef iterator::reference                  reference;
    00589     typedef iterator::index_reference            index_reference;
    00590     typedef iterator::pointer                    pointer;
    00591     typedef iterator::difference_type            difference_type;
    00592     typedef iterator::row_iterator               row_iterator;
    00593     typedef iterator::column_iterator            column_iterator;
    00594     typedef VectorAccessor<FFTWComplex>          default_accessor;
    00595     typedef VectorAccessor<FFTWComplex>          DefaultAccessor;
    00596 };
    00597
    00598     /** Complex (FFTWComplex) image.
    00599         It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
    00600         their const counterparts to access the data.
    00601
    00602         <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00603         <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00604         Namespace: vigra
    00605     */
    00606 typedef BasicImage<FFTWComplex> FFTWComplexImage;
    00607
    00608 //@}
    00609
    00610 /********************************************************/
    00611 /*                                                      */
    00612 /*                  FFTWComplex-Accessors               */
    00613 /*                                                      */
    00614 /********************************************************/
    00615
    00616 /** \addtogroup DataAccessors
    00617 */
    00618 //@{
    00619 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex
    00620
    00621     Encapsulate access to pixels of type FFTWComplex
    00622 */
    00623 //@{
    00624     /** Encapsulate access to the the real part of a complex number.
    00625
    00626     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00627     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00628     Namespace: vigra
    00629     */
    00630 class FFTWRealAccessor
    00631 {
    00632   public:
    00633
    00634         /// The accessor's value type.
    00635     typedef fftw_real value_type;
    00636
    00637         /// Read real part at iterator position.
    00638     template <class ITERATOR>
    00639     value_type operator()(ITERATOR const & i) const {
    00640         return (*i).re();
    00641     }
    00642
    00643         /// Read real part at offset from iterator position.
    00644     template <class ITERATOR, class DIFFERENCE>
    00645     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00646         return i[d].re();
    00647     }
    00648
    00649         /// Write real part at iterator position.
    00650     template <class ITERATOR>
    00651     void set(value_type const & v, ITERATOR const & i) const {
    00652         (*i).re()= v;
    00653     }
    00654
    00655         /// Write real part at offset from iterator position.
    00656     template <class ITERATOR, class DIFFERENCE>
    00657     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
    00658         i[d].re()= v;
    00659     }
    00660 };
    00661
    00662     /** Encapsulate access to the the imaginary part of a complex number.
    00663
    00664     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00665     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00666     Namespace: vigra
    00667     */
    00668 class FFTWImaginaryAccessor
    00669 {
    00670   public:
    00671         /// The accessor's value type.
    00672     typedef fftw_real value_type;
    00673
    00674         /// Read imaginary part at iterator position.
    00675     template <class ITERATOR>
    00676     value_type operator()(ITERATOR const & i) const {
    00677         return (*i).im();
    00678     }
    00679
    00680         /// Read imaginary part at offset from iterator position.
    00681     template <class ITERATOR, class DIFFERENCE>
    00682     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00683         return i[d].im();
    00684     }
    00685
    00686         /// Write imaginary part at iterator position.
    00687     template <class ITERATOR>
    00688     void set(value_type const & v, ITERATOR const & i) const {
    00689         (*i).im()= v;
    00690     }
    00691
    00692         /// Write imaginary part at offset from iterator position.
    00693     template <class ITERATOR, class DIFFERENCE>
    00694     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
    00695         i[d].im()= v;
    00696     }
    00697 };
    00698
    00699     /** Write a real number into a complex one. The imaginary part is set
    00700         to 0.
    00701
    00702     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00703     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00704     Namespace: vigra
    00705     */
    00706 class FFTWWriteRealAccessor: public FFTWRealAccessor
    00707 {
    00708   public:
    00709         /// The accessor's value type.
    00710     typedef fftw_real value_type;
    00711
    00712         /** Write real number at iterator position. Set imaginary part
    00713             to 0.
    00714         */
    00715     template <class ITERATOR>
    00716     void set(value_type const & v, ITERATOR const & i) const {
    00717         (*i).re()= v;
    00718         (*i).im()= 0;
    00719     }
    00720
    00721         /** Write real number at offset from iterator position. Set imaginary part
    00722             to 0.
    00723         */
    00724     template <class ITERATOR, class DIFFERENCE>
    00725     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
    00726         i[d].re()= v;
    00727         i[d].im()= 0;
    00728     }
    00729 };
    00730
    00731     /** Calculate magnitude of complex number on the fly.
    00732
    00733     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00734     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00735     Namespace: vigra
    00736     */
    00737 class FFTWMagnitudeAccessor
    00738 {
    00739   public:
    00740         /// The accessor's value type.
    00741     typedef fftw_real value_type;
    00742
    00743         /// Read magnitude at iterator position.
    00744     template <class ITERATOR>
    00745     value_type operator()(ITERATOR const & i) const {
    00746         return (*i).magnitude();
    00747     }
    00748
    00749         /// Read magnitude at offset from iterator position.
    00750     template <class ITERATOR, class DIFFERENCE>
    00751     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00752         return (i[d]).magnitude();
    00753     }
    00754 };
    00755
    00756     /** Calculate phase of complex number on the fly.
    00757
    00758     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>" (for FFTW 3) or<br>
    00759     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>" (for deprecated FFTW 2)<br>
    00760     Namespace: vigra
    00761     */
    00762 class FFTWPhaseAccessor
    00763 {
    00764   public:
    00765         /// The accessor's value type.
    00766     typedef fftw_real value_type;
    00767
    00768         /// Read phase at iterator position.
    00769     template <class ITERATOR>
    00770     value_type operator()(ITERATOR const & i) const {
    00771         return (*i).phase();
    00772     }
    00773
    00774         /// Read phase at offset from iterator position.
    00775     template <class ITERATOR, class DIFFERENCE>
    00776     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00777         return (i[d]).phase();
    00778     }
    00779 };
    00780
    00781 //@}
    00782 //@}
    00783
    00784 /********************************************************/
    00785 /*                                                      */
    00786 /*                    Fourier Transform                 */
    00787 /*                                                      */
    00788 /********************************************************/
    00789
    00790 /** \addtogroup FourierTransform Fast Fourier Transform
    00791
    00792     This documentation describes the VIGRA interface to FFTW 3. There is also an
    00793     \link FourierTransformFFTW2 interface to the older version FFTW 2\endlink
    00794
    00795     VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
    00796     Transform</a> package to perform Fourier transformations. VIGRA
    00797     provides a wrapper for FFTW's complex number type (FFTWComplex),
    00798     but FFTW's functions are used verbatim. If the image is stored as
    00799     a FFTWComplexImage, the simplest call to an FFT function is like this:
    00800
    00801     \code
    00802     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
    00803     ... // fill image with data
    00804
    00805     // create a plan with estimated performance optimization
    00806     fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
    00807                                 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
    00808                                 FFTW_FORWARD, FFTW_ESTIMATE );
    00809     // calculate FFT (this can be repeated as often as needed,
    00810     //                with fresh data written into the source array)
    00811     fftw_execute(forwardPlan);
    00812
    00813     // release the plan memory
    00814     fftw_destroy_plan(forwardPlan);
    00815
    00816     // likewise for the inverse transform
    00817     fftw_plan backwardPlan = fftw_plan_dft_2d(height, width,
    00818                                  (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(),
    00819                                  FFTW_BACKWARD, FFTW_ESTIMATE);
    00820     fftw_execute(backwardPlan);
    00821     fftw_destroy_plan(backwardPlan);
    00822
    00823     // do not forget to normalize the result according to the image size
    00824     transformImage(srcImageRange(spatial), destImage(spatial),
    00825                    std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height));
    00826     \endcode
    00827
    00828     Note that in the creation of a plan, the height must be given
    00829     first. Note also that <TT>spatial.begin()</TT> may only be passed
    00830     to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the
    00831     entire image. When you want to restrict operation to an ROI, you
    00832     can create a copy of the ROI in an image of appropriate size, or
    00833     you may use the Guru interface to FFTW.
    00834
    00835     More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
    00836
    00837     FFTW produces fourier images that have the DC component (the
    00838     origin of the Fourier space) in the upper left corner. Often, one
    00839     wants the origin in the center of the image, so that frequencies
    00840     always increase towards the border of the image. This can be
    00841     achieved by calling \ref moveDCToCenter(). The inverse
    00842     transformation is done by \ref moveDCToUpperLeft().
    00843
    00844     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
    00845     Namespace: vigra
    00846 */
    00847
    00848 /** \addtogroup FourierTransform
    00849 */
    00850 //@{
    00851
    00852 /********************************************************/
    00853 /*                                                      */
    00854 /*                     moveDCToCenter                   */
    00855 /*                                                      */
    00856 /********************************************************/
    00857
    00858 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
    00859           in the image center.
    00860
    00861     FFTW produces fourier images where the DC component (origin of
    00862     fourier space) is located in the upper left corner of the
    00863     image. The quadrants are placed like this (using a 4x4 image for
    00864     example):
    00865
    00866     \code
    00867             DC 4 3 3
    00868              4 4 3 3
    00869              1 1 2 2
    00870              1 1 2 2
    00871     \endcode
    00872
    00873     After applying the function, the quadrants are at their usual places:
    00874
    00875     \code
    00876             2 2  1 1
    00877             2 2  1 1
    00878             3 3 DC 4
    00879             3 3  4 4
    00880     \endcode
    00881
    00882     This transformation can be reversed by \ref moveDCToUpperLeft().
    00883     Note that the transformation must not be executed in place - input
    00884     and output images must be different.
    00885
    00886     <b> Declarations:</b>
    00887
    00888     pass arguments explicitly:
    00889     \code
    00890     namespace vigra {
    00891         template <class SrcImageIterator, class SrcAccessor,
    00892           class DestImageIterator, class DestAccessor>
    00893         void moveDCToCenter(SrcImageIterator src_upperleft,
    00894                                SrcImageIterator src_lowerright, SrcAccessor sa,
    00895                                DestImageIterator dest_upperleft, DestAccessor da);
    00896     }
    00897     \endcode
    00898
    00899
    00900     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00901     \code
    00902     namespace vigra {
    00903         template <class SrcImageIterator, class SrcAccessor,
    00904                   class DestImageIterator, class DestAccessor>
    00905         void moveDCToCenter(
    00906             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00907             pair<DestImageIterator, DestAccessor> dest);
    00908     }
    00909     \endcode
    00910
    00911     <b> Usage:</b>
    00912
    00913         <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
    00914         Namespace: vigra
    00915
    00916     \code
    00917     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
    00918     ... // fill image with data
    00919
    00920     // create a plan with estimated performance optimization
    00921     fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
    00922                                 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
    00923                                 FFTW_FORWARD, FFTW_ESTIMATE );
    00924     // calculate FFT
    00925     fftw_execute(forwardPlan);
    00926
    00927     vigra::FFTWComplexImage rearrangedFourier(width, height);
    00928     moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
    00929
    00930     // delete the plan
    00931     fftw_destroy_plan(forwardPlan);
    00932     \endcode
    00933 */
    00934 template <class SrcImageIterator, class SrcAccessor,
    00935           class DestImageIterator, class DestAccessor>
    00936 void moveDCToCenter(SrcImageIterator src_upperleft,
    00937                                SrcImageIterator src_lowerright, SrcAccessor sa,
    00938                                DestImageIterator dest_upperleft, DestAccessor da)
    00939 {
    00940     int w= src_lowerright.x - src_upperleft.x;
    00941     int h= src_lowerright.y - src_upperleft.y;
    00942     int w1 = w/2;
    00943     int h1 = h/2;
    00944     int w2 = (w+1)/2;
    00945     int h2 = (h+1)/2;
    00946
    00947     // 2. Quadrant  zum 4.
    00948     copyImage(srcIterRange(src_upperleft,
    00949                            src_upperleft  + Diff2D(w2, h2), sa),
    00950               destIter    (dest_upperleft + Diff2D(w1, h1), da));
    00951
    00952     // 4. Quadrant zum 2.
    00953     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
    00954                            src_lowerright, sa),
    00955               destIter    (dest_upperleft, da));
    00956
    00957     // 1. Quadrant zum 3.
    00958     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
    00959                            src_upperleft  + Diff2D(w,  h2), sa),
    00960               destIter    (dest_upperleft + Diff2D(0,  h1), da));
    00961
    00962     // 3. Quadrant zum 1.
    00963     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
    00964                            src_upperleft  + Diff2D(w2, h), sa),
    00965               destIter    (dest_upperleft + Diff2D(w1, 0), da));
    00966 }
    00967
    00968 template <class SrcImageIterator, class SrcAccessor,
    00969           class DestImageIterator, class DestAccessor>
    00970 inline void moveDCToCenter(
    00971     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00972     pair<DestImageIterator, DestAccessor> dest)
    00973 {
    00974     moveDCToCenter(src.first, src.second, src.third,
    00975                           dest.first, dest.second);
    00976 }
    00977
    00978 /********************************************************/
    00979 /*                                                      */
    00980 /*                   moveDCToUpperLeft                  */
    00981 /*                                                      */
    00982 /********************************************************/
    00983
    00984 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
    00985           in the image's upper left.
    00986
    00987      This function is the inversion of \ref moveDCToCenter(). See there
    00988      for declarations and a usage example.
    00989
    00990     <b> Declarations:</b>
    00991
    00992     pass arguments explicitly:
    00993     \code
    00994     namespace vigra {
    00995         template <class SrcImageIterator, class SrcAccessor,
    00996           class DestImageIterator, class DestAccessor>
    00997         void moveDCToUpperLeft(SrcImageIterator src_upperleft,
    00998                                SrcImageIterator src_lowerright, SrcAccessor sa,
    00999                                DestImageIterator dest_upperleft, DestAccessor da);
    01000     }
    01001     \endcode
    01002
    01003
    01004     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01005     \code
    01006     namespace vigra {
    01007         template <class SrcImageIterator, class SrcAccessor,
    01008                   class DestImageIterator, class DestAccessor>
    01009         void moveDCToUpperLeft(
    01010             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01011             pair<DestImageIterator, DestAccessor> dest);
    01012     }
    01013     \endcode
    01014 */
    01015 template <class SrcImageIterator, class SrcAccessor,
    01016           class DestImageIterator, class DestAccessor>
    01017 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
    01018                                SrcImageIterator src_lowerright, SrcAccessor sa,
    01019                                DestImageIterator dest_upperleft, DestAccessor da)
    01020 {
    01021     int w= src_lowerright.x - src_upperleft.x;
    01022     int h= src_lowerright.y - src_upperleft.y;
    01023     int w2 = w/2;
    01024     int h2 = h/2;
    01025     int w1 = (w+1)/2;
    01026     int h1 = (h+1)/2;
    01027
    01028     // 2. Quadrant  zum 4.
    01029     copyImage(srcIterRange(src_upperleft,
    01030                            src_upperleft  + Diff2D(w2, h2), sa),
    01031               destIter    (dest_upperleft + Diff2D(w1, h1), da));
    01032
    01033     // 4. Quadrant zum 2.
    01034     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
    01035                            src_lowerright, sa),
    01036               destIter    (dest_upperleft, da));
    01037
    01038     // 1. Quadrant zum 3.
    01039     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
    01040                            src_upperleft  + Diff2D(w,  h2), sa),
    01041               destIter    (dest_upperleft + Diff2D(0,  h1), da));
    01042
    01043     // 3. Quadrant zum 1.
    01044     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
    01045                            src_upperleft  + Diff2D(w2, h), sa),
    01046               destIter    (dest_upperleft + Diff2D(w1, 0), da));
    01047 }
    01048
    01049 template <class SrcImageIterator, class SrcAccessor,
    01050           class DestImageIterator, class DestAccessor>
    01051 inline void moveDCToUpperLeft(
    01052     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01053     pair<DestImageIterator, DestAccessor> dest)
    01054 {
    01055     moveDCToUpperLeft(src.first, src.second, src.third,
    01056                                           dest.first, dest.second);
    01057 }
    01058
    01059 /********************************************************/
    01060 /*                                                      */
    01061 /*                   applyFourierFilter                 */
    01062 /*                                                      */
    01063 /********************************************************/
    01064
    01065 /** \brief Apply a filter (defined in the frequency domain) to an image.
    01066
    01067     After transferring the image into the frequency domain, it is
    01068     multiplied pixel-wise with the filter and transformed back. The
    01069     result is always put into the given FFTWComplexImage
    01070     <TT>destImg</TT> which must have the right size. Finally, the
    01071     result will be normalized to compensate for the two FFTs.
    01072
    01073     The input and filter images can be scalar images (more precisely,
    01074     <TT>SrcAccessor::value_type</TT> must be scalar) or
    01075     FFTWComplexImages (in this case, one conversion is saved).
    01076
    01077     The DC entry of the filter must be in the upper left, which is the
    01078     position where FFTW expects it (see \ref moveDCToUpperLeft()).
    01079
    01080     You can optionally pass two fftwnd_plans as last parameters, the
    01081     forward plan and the in-place backward plan. Both must have been
    01082     created for the right image size (see sample code).
    01083
    01084     <b> Declarations:</b>
    01085
    01086     pass arguments explicitly:
    01087     \code
    01088     namespace vigra {
    01089         template <class SrcImageIterator, class SrcAccessor,
    01090             class FilterImageIterator, class FilterAccessor>
    01091         void applyFourierFilter(SrcImageIterator srcUpperLeft,
    01092                                 SrcImageIterator srcLowerRight, SrcAccessor sa,
    01093                                 FilterImageIterator filterUpperLeft, FilterAccessor fa,
    01094                                 FFTWComplexImage & destImg)
    01095
    01096         template <class SrcImageIterator, class SrcAccessor,
    01097             class FilterImageIterator, class FilterAccessor>
    01098         void applyFourierFilter(SrcImageIterator srcUpperLeft,
    01099                                 SrcImageIterator srcLowerRight, SrcAccessor sa,
    01100                                 FilterImageIterator filterUpperLeft, FilterAccessor fa,
    01101                                 FFTWComplexImage & destImg,
    01102                                 const fftwnd_plan & forwardPlan, const fftwnd_plan & backwardPlan)
    01103     }
    01104     \endcode
    01105
    01106     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01107     \code
    01108     namespace vigra {
    01109         template <class SrcImageIterator, class SrcAccessor,
    01110             class FilterImageIterator, class FilterAccessor>
    01111         inline
    01112         void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01113                                 pair<FilterImageIterator, FilterAccessor> filter,
    01114                                 FFTWComplexImage &destImg)
    01115
    01116         template <class SrcImageIterator, class SrcAccessor,
    01117             class FilterImageIterator, class FilterAccessor>
    01118         inline
    01119         void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01120                                 pair<FilterImageIterator, FilterAccessor> filter,
    01121                                 FFTWComplexImage &destImg,
    01122                                 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    01123     }
    01124     \endcode
    01125
    01126     <b> Usage:</b>
    01127
    01128     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
    01129     Namespace: vigra
    01130
    01131     \code
    01132     // create a Gaussian filter in Fourier space
    01133     vigra::FImage gaussFilter(w, h), filter(w, h);
    01134     for(int y=0; y<h; ++y)
    01135         for(int x=0; x<w; ++x)
    01136         {
    01137             xx = float(x - w / 2) / w;
    01138             yy = float(y - h / 2) / h;
    01139
    01140             gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
    01141         }
    01142
    01143     // applyFourierFilter() expects the filter's DC in the upper left
    01144     moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
    01145
    01146     vigra::FFTWComplexImage result(w, h);
    01147
    01148     vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
    01149     \endcode
    01150
    01151     For inspection of the result, \ref FFTWMagnitudeAccessor might be
    01152     useful. If you want to apply the same filter repeatedly, it may be more
    01153     efficient to use the FFTW functions directly with FFTW plans optimized
    01154     for good performance.
    01155 */
    01156 template <class SrcImageIterator, class SrcAccessor,
    01157           class FilterImageIterator, class FilterAccessor,
    01158           class DestImageIterator, class DestAccessor>
    01159 inline
    01160 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01161                         pair<FilterImageIterator, FilterAccessor> filter,
    01162                         pair<DestImageIterator, DestAccessor> dest)
    01163 {
    01164     applyFourierFilter(src.first, src.second, src.third,
    01165                        filter.first, filter.second,
    01166                        dest.first, dest.second);
    01167 }
    01168
    01169 template <class SrcImageIterator, class SrcAccessor,
    01170           class FilterImageIterator, class FilterAccessor,
    01171           class DestImageIterator, class DestAccessor>
    01172 void applyFourierFilter(SrcImageIterator srcUpperLeft,
    01173                         SrcImageIterator srcLowerRight, SrcAccessor sa,
    01174                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
    01175                         DestImageIterator destUpperLeft, DestAccessor da)
    01176 {
    01177     // copy real input images into a complex one...
    01178     int w= srcLowerRight.x - srcUpperLeft.x;
    01179     int h= srcLowerRight.y - srcUpperLeft.y;
    01180
    01181     FFTWComplexImage workImage(w, h);
    01182     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    01183               destImage(workImage, FFTWWriteRealAccessor()));
    01184
    01185     // ...and call the impl
    01186     FFTWComplexImage const & cworkImage = workImage;
    01187     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    01188                            filterUpperLeft, fa,
    01189                            destUpperLeft, da);
    01190 }
    01191
    01192 template <class FilterImageIterator, class FilterAccessor,
    01193           class DestImageIterator, class DestAccessor>
    01194 inline
    01195 void applyFourierFilter(
    01196     FFTWComplexImage::const_traverser srcUpperLeft,
    01197     FFTWComplexImage::const_traverser srcLowerRight,
    01198     FFTWComplexImage::ConstAccessor sa,
    01199     FilterImageIterator filterUpperLeft, FilterAccessor fa,
    01200     DestImageIterator destUpperLeft, DestAccessor da)
    01201 {
    01202     int w = srcLowerRight.x - srcUpperLeft.x;
    01203     int h = srcLowerRight.y - srcUpperLeft.y;
    01204
    01205     // test for right memory layout (fftw expects a 2*width*height floats array)
    01206     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
    01207         applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
    01208                                filterUpperLeft, fa,
    01209                                destUpperLeft, da);
    01210     else
    01211     {
    01212         FFTWComplexImage workImage(w, h);
    01213         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    01214                   destImage(workImage));
    01215
    01216         FFTWComplexImage const & cworkImage = workImage;
    01217         applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    01218                                filterUpperLeft, fa,
    01219                                destUpperLeft, da);
    01220     }
    01221 }
    01222
    01223 template <class FilterImageIterator, class FilterAccessor,
    01224           class DestImageIterator, class DestAccessor>
    01225 void applyFourierFilterImpl(
    01226     FFTWComplexImage::const_traverser srcUpperLeft,
    01227     FFTWComplexImage::const_traverser srcLowerRight,
    01228     FFTWComplexImage::ConstAccessor sa,
    01229     FilterImageIterator filterUpperLeft, FilterAccessor fa,
    01230     DestImageIterator destUpperLeft, DestAccessor da)
    01231 {
    01232     int w = srcLowerRight.x - srcUpperLeft.x;
    01233     int h = srcLowerRight.y - srcUpperLeft.y;
    01234
    01235     FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
    01236
    01237     // FFT from srcImage to complexResultImg
    01238     fftw_plan forwardPlan=
    01239         fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
    01240                                (fftw_complex *)complexResultImg.begin(),
    01241                                FFTW_FORWARD, FFTW_ESTIMATE );
    01242     fftw_execute(forwardPlan);
    01243     fftw_destroy_plan(forwardPlan);
    01244
    01245     // convolve in freq. domain (in complexResultImg)
    01246     combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
    01247                      destImage(complexResultImg), std::multiplies<FFTWComplex>());
    01248
    01249     // FFT back into spatial domain (inplace in complexResultImg)
    01250     fftw_plan backwardPlan=
    01251         fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
    01252                                (fftw_complex *)complexResultImg.begin(),
    01253                                FFTW_BACKWARD, FFTW_ESTIMATE);
    01254     fftw_execute(backwardPlan);
    01255     fftw_destroy_plan(backwardPlan);
    01256
    01257     typedef typename
    01258         NumericTraits<typename DestAccessor::value_type>::isScalar
    01259         isScalarResult;
    01260
    01261     // normalization (after FFTs), maybe stripping imaginary part
    01262     applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
    01263                                         isScalarResult());
    01264 }
    01265
    01266 template <class DestImageIterator, class DestAccessor>
    01267 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
    01268                                          DestImageIterator destUpperLeft,
    01269                                          DestAccessor da,
    01270                                          VigraFalseType)
    01271 {
    01272     double normFactor= 1.0/(srcImage.width() * srcImage.height());
    01273
    01274     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
    01275     {
    01276         DestImageIterator dIt= destUpperLeft;
    01277         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
    01278         {
    01279             da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
    01280             da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
    01281         }
    01282     }
    01283 }
    01284
    01285 inline
    01286 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
    01287         FFTWComplexImage::traverser destUpperLeft,
    01288         FFTWComplexImage::Accessor da,
    01289         VigraFalseType)
    01290 {
    01291     transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
    01292                    linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height())));
    01293 }
    01294
    01295 template <class DestImageIterator, class DestAccessor>
    01296 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
    01297                                          DestImageIterator destUpperLeft,
    01298                                          DestAccessor da,
    01299                                          VigraTrueType)
    01300 {
    01301     double normFactor= 1.0/(srcImage.width() * srcImage.height());
    01302
    01303     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
    01304     {
    01305         DestImageIterator dIt= destUpperLeft;
    01306         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
    01307             da.set(srcImage(x, y).re()*normFactor, dIt);
    01308     }
    01309 }
    01310
    01311 /**********************************************************/
    01312 /*                                                        */
    01313 /*                applyFourierFilterFamily                */
    01314 /*                                                        */
    01315 /**********************************************************/
    01316
    01317 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
    01318
    01319     This provides the same functionality as \ref applyFourierFilter(),
    01320     but applying several filters at once allows to avoid
    01321     repeated Fourier transforms of the source image.
    01322
    01323     Filters and result images must be stored in \ref vigra::ImageArray data
    01324     structures. In contrast to \ref applyFourierFilter(), this function adjusts
    01325     the size of the result images and the the length of the array.
    01326
    01327     <b> Declarations:</b>
    01328
    01329     pass arguments explicitly:
    01330     \code
    01331     namespace vigra {
    01332         template <class SrcImageIterator, class SrcAccessor, class FilterType>
    01333         void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
    01334                                       SrcImageIterator srcLowerRight, SrcAccessor sa,
    01335                                       const ImageArray<FilterType> &filters,
    01336                                       ImageArray<FFTWComplexImage> &results)
    01337     }
    01338     \endcode
    01339
    01340     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01341     \code
    01342     namespace vigra {
    01343         template <class SrcImageIterator, class SrcAccessor, class FilterType>
    01344         inline
    01345         void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01346                                       const ImageArray<FilterType> &filters,
    01347                                       ImageArray<FFTWComplexImage> &results)
    01348     }
    01349     \endcode
    01350
    01351     <b> Usage:</b>
    01352
    01353     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
    01354     Namespace: vigra
    01355
    01356     \code
    01357     // assuming the presence of a real-valued image named "image" to
    01358     // be filtered in this example
    01359
    01360     vigra::ImageArray<vigra::FImage> filters(16, image.size());
    01361
    01362     for (int i=0; i<filters.size(); i++)
    01363          // create some meaningful filters here
    01364          createMyFilterOfScale(i, destImage(filters[i]));
    01365
    01366     vigra::ImageArray<vigra::FFTWComplexImage> results();
    01367
    01368     vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
    01369     \endcode
    01370 */
    01371 template <class SrcImageIterator, class SrcAccessor,
    01372           class FilterType, class DestImage>
    01373 inline
    01374 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01375                               const ImageArray<FilterType> &filters,
    01376                               ImageArray<DestImage> &results)
    01377 {
    01378     applyFourierFilterFamily(src.first, src.second, src.third,
    01379                              filters, results);
    01380 }
    01381
    01382 template <class SrcImageIterator, class SrcAccessor,
    01383           class FilterType, class DestImage>
    01384 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
    01385                               SrcImageIterator srcLowerRight, SrcAccessor sa,
    01386                               const ImageArray<FilterType> &filters,
    01387                               ImageArray<DestImage> &results)
    01388 {
    01389     int w= srcLowerRight.x - srcUpperLeft.x;
    01390     int h= srcLowerRight.y - srcUpperLeft.y;
    01391
    01392     FFTWComplexImage workImage(w, h);
    01393     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    01394               destImage(workImage, FFTWWriteRealAccessor()));
    01395
    01396     FFTWComplexImage const & cworkImage = workImage;
    01397     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    01398                                  filters, results);
    01399 }
    01400
    01401 template <class FilterType, class DestImage>
    01402 inline
    01403 void applyFourierFilterFamily(
    01404     FFTWComplexImage::const_traverser srcUpperLeft,
    01405     FFTWComplexImage::const_traverser srcLowerRight,
    01406     FFTWComplexImage::ConstAccessor sa,
    01407     const ImageArray<FilterType> &filters,
    01408     ImageArray<DestImage> &results)
    01409 {
    01410     int w= srcLowerRight.x - srcUpperLeft.x;
    01411
    01412     // test for right memory layout (fftw expects a 2*width*height floats array)
    01413     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
    01414         applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
    01415                                      filters, results);
    01416     else
    01417     {
    01418         int h = srcLowerRight.y - srcUpperLeft.y;
    01419         FFTWComplexImage workImage(w, h);
    01420         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    01421                   destImage(workImage));
    01422
    01423         FFTWComplexImage const & cworkImage = workImage;
    01424         applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    01425                                      filters, results);
    01426     }
    01427 }
    01428
    01429 template <class FilterType, class DestImage>
    01430 void applyFourierFilterFamilyImpl(
    01431     FFTWComplexImage::const_traverser srcUpperLeft,
    01432     FFTWComplexImage::const_traverser srcLowerRight,
    01433     FFTWComplexImage::ConstAccessor sa,
    01434     const ImageArray<FilterType> &filters,
    01435     ImageArray<DestImage> &results)
    01436 {
    01437     // make sure the filter images have the right dimensions
    01438     vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
    01439                        "applyFourierFilterFamily called with src image size != filters.imageSize()!");
    01440
    01441     // make sure the result image array has the right dimensions
    01442     results.resize(filters.size());
    01443     results.resizeImages(filters.imageSize());
    01444
    01445     // FFT from srcImage to freqImage
    01446     int w= srcLowerRight.x - srcUpperLeft.x;
    01447     int h= srcLowerRight.y - srcUpperLeft.y;
    01448
    01449     FFTWComplexImage freqImage(w, h);
    01450     FFTWComplexImage result(w, h);
    01451
    01452     fftw_plan forwardPlan=
    01453         fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
    01454                                (fftw_complex *)freqImage.begin(),
    01455                                FFTW_FORWARD, FFTW_ESTIMATE );
    01456     fftw_execute(forwardPlan);
    01457     fftw_destroy_plan(forwardPlan);
    01458
    01459     fftw_plan backwardPlan=
    01460         fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
    01461                                (fftw_complex *)result.begin(),
    01462                                FFTW_BACKWARD, FFTW_ESTIMATE );
    01463     typedef typename
    01464         NumericTraits<typename DestImage::Accessor::value_type>::isScalar
    01465         isScalarResult;
    01466
    01467     // convolve with filters in freq. domain
    01468     for (unsigned int i= 0;  i < filters.size(); i++)
    01469     {
    01470         combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
    01471                          destImage(result), std::multiplies<FFTWComplex>());
    01472
    01473         // FFT back into spatial domain (inplace in destImage)
    01474         fftw_execute(backwardPlan);
    01475
    01476         // normalization (after FFTs), maybe stripping imaginary part
    01477         applyFourierFilterImplNormalization(result,
    01478                                             results[i].upperLeft(), results[i].accessor(),
    01479                                             isScalarResult());
    01480     }
    01481     fftw_destroy_plan(backwardPlan);
    01482 }
    01483
    01484 /********************************************************/
    01485 /*                                                      */
    01486 /*                fourierTransformReal                  */
    01487 /*                                                      */
    01488 /********************************************************/
    01489
    01490 /** \brief Real Fourier transforms for even and odd boundary conditions
    01491            (aka. cosine and sine transforms).
    01492
    01493
    01494     If the image is real and has even symmetry, its Fourier transform
    01495     is also real and has even symmetry. The Fourier transform of a real image with odd
    01496     symmetry is imaginary and has odd symmetry. In either case, only about a quarter
    01497     of the pixels need to be stored because the rest can be calculated from the symmetry
    01498     properties. This is especially useful, if the original image is implicitly assumed
    01499     to have reflective or anti-reflective boundary conditions. Then the "negative"
    01500     pixel locations are defined as
    01501
    01502     \code
    01503     even (reflective boundary conditions):      f[-x] = f[x]     (x = 1,...,N-1)
    01504     odd (anti-reflective boundary conditions):  f[-1] = 0
    01505                                                 f[-x] = -f[x-2]  (x = 2,...,N-1)
    01506     \endcode
    01507
    01508     end similar at the other boundary (see the FFTW documentation for details).
    01509     This has the advantage that more efficient Fourier transforms that use only
    01510     real numbers can be implemented. These are also known as cosine and sine transforms
    01511     respectively.
    01512
    01513     If you use the odd transform it is important to note that in the Fourier domain,
    01514     the DC component is always zero and is therefore dropped from the data structure.
    01515     This means that index 0 in an odd symmetric Fourier domain image refers to
    01516     the <i>first</i> harmonic. This is especially important if an image is first
    01517     cosine transformed (even symmetry), then in the Fourier domain multiplied
    01518     with an odd symmetric filter (e.g. a first derivative) and finally transformed
    01519     back to the spatial domain with a sine transform (odd symmetric). For this to work
    01520     properly the image must be shifted left or up by one pixel (depending on whether
    01521     the x- or y-axis is odd symmetric) before the inverse transform can be applied.
    01522     (see example below).
    01523
    01524     The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
    01525     where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
    01526     whether the x- and y-axis is to be transformed using even or odd symmetry.
    01527     The same functions can be used for both the forward and inverse transforms,
    01528     only the normalization changes. For signal processing, the following
    01529     normalization factors are most appropriate:
    01530
    01531     \code
    01532                           forward             inverse
    01533     ------------------------------------------------------------
    01534     X even, Y even           1.0         4.0 * (w-1) * (h-1)
    01535     X even, Y odd           -1.0        -4.0 * (w-1) * (h+1)
    01536     X odd,  Y even          -1.0        -4.0 * (w+1) * (h-1)
    01537     X odd,  Y odd            1.0         4.0 * (w+1) * (h+1)
    01538     \endcode
    01539
    01540     where <tt>w</tt> and <tt>h</tt> denote the image width and height.
    01541
    01542     <b> Declarations:</b>
    01543
    01544     pass arguments explicitly:
    01545     \code
    01546     namespace vigra {
    01547         template <class SrcTraverser, class SrcAccessor,
    01548                   class DestTraverser, class DestAccessor>
    01549         void
    01550         fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
    01551                                DestTraverser dul, DestAccessor dest, fftw_real norm);
    01552
    01553         fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
    01554     }
    01555     \endcode
    01556
    01557
    01558     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01559     \code
    01560     namespace vigra {
    01561         template <class SrcTraverser, class SrcAccessor,
    01562                   class DestTraverser, class DestAccessor>
    01563         void
    01564         fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
    01565                                pair<DestTraverser, DestAccessor> dest, fftw_real norm);
    01566
    01567         fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
    01568     }
    01569     \endcode
    01570
    01571     <b> Usage:</b>
    01572
    01573         <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
    01574         Namespace: vigra
    01575
    01576     \code
    01577     vigra::FImage spatial(width,height), fourier(width,height);
    01578     ... // fill image with data
    01579
    01580     // forward cosine transform == reflective boundary conditions
    01581     fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
    01582
    01583     // multiply with a first derivative of Gaussian in x-direction
    01584     for(int y = 0; y < height; ++y)
    01585     {
    01586         for(int x = 1; x < width; ++x)
    01587         {
    01588             double dx = x * M_PI / (width - 1);
    01589             double dy = y * M_PI / (height - 1);
    01590             fourier(x-1, y) = fourier(x, y) * dx * exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
    01591         }
    01592         fourier(width-1, y) = 0.0;
    01593     }
    01594
    01595     // inverse transform -- odd symmetry in x-direction, even in y,
    01596     //                      due to symmetry of the filter
    01597     fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
    01598                            (fftw_real)-4.0 * (width+1) * (height-1));
    01599     \endcode
    01600 */
    01601 template <class SrcTraverser, class SrcAccessor,
    01602           class DestTraverser, class DestAccessor>
    01603 inline void
    01604 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
    01605                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
    01606 {
    01607     fourierTransformRealEE(src.first, src.second, src.third,
    01608                                    dest.first, dest.second, norm);
    01609 }
    01610
    01611 template <class SrcTraverser, class SrcAccessor,
    01612           class DestTraverser, class DestAccessor>
    01613 inline void
    01614 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
    01615                                DestTraverser dul, DestAccessor dest, fftw_real norm)
    01616 {
    01617     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01618                                       norm, FFTW_REDFT00, FFTW_REDFT00);
    01619 }
    01620
    01621 template <class DestTraverser, class DestAccessor>
    01622 inline void
    01623 fourierTransformRealEE(
    01624          FFTWRealImage::const_traverser sul,
    01625          FFTWRealImage::const_traverser slr,
    01626          FFTWRealImage::Accessor src,
    01627          DestTraverser dul, DestAccessor dest, fftw_real norm)
    01628 {
    01629     int w = slr.x - sul.x;
    01630
    01631     // test for right memory layout (fftw expects a width*height fftw_real array)
    01632     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
    01633         fourierTransformRealImpl(sul, slr, dul, dest,
    01634                                  norm, FFTW_REDFT00, FFTW_REDFT00);
    01635     else
    01636         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01637                                  norm, FFTW_REDFT00, FFTW_REDFT00);
    01638 }
    01639
    01640 /********************************************************************/
    01641
    01642 template <class SrcTraverser, class SrcAccessor,
    01643           class DestTraverser, class DestAccessor>
    01644 inline void
    01645 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
    01646                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
    01647 {
    01648     fourierTransformRealOE(src.first, src.second, src.third,
    01649                                    dest.first, dest.second, norm);
    01650 }
    01651
    01652 template <class SrcTraverser, class SrcAccessor,
    01653           class DestTraverser, class DestAccessor>
    01654 inline void
    01655 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
    01656                                DestTraverser dul, DestAccessor dest, fftw_real norm)
    01657 {
    01658     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01659                                       norm, FFTW_RODFT00, FFTW_REDFT00);
    01660 }
    01661
    01662 template <class DestTraverser, class DestAccessor>
    01663 inline void
    01664 fourierTransformRealOE(
    01665          FFTWRealImage::const_traverser sul,
    01666          FFTWRealImage::const_traverser slr,
    01667          FFTWRealImage::Accessor src,
    01668          DestTraverser dul, DestAccessor dest, fftw_real norm)
    01669 {
    01670     int w = slr.x - sul.x;
    01671
    01672     // test for right memory layout (fftw expects a width*height fftw_real array)
    01673     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
    01674         fourierTransformRealImpl(sul, slr, dul, dest,
    01675                                  norm, FFTW_RODFT00, FFTW_REDFT00);
    01676     else
    01677         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01678                                  norm, FFTW_RODFT00, FFTW_REDFT00);
    01679 }
    01680
    01681 /********************************************************************/
    01682
    01683 template <class SrcTraverser, class SrcAccessor,
    01684           class DestTraverser, class DestAccessor>
    01685 inline void
    01686 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
    01687                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
    01688 {
    01689     fourierTransformRealEO(src.first, src.second, src.third,
    01690                                    dest.first, dest.second, norm);
    01691 }
    01692
    01693 template <class SrcTraverser, class SrcAccessor,
    01694           class DestTraverser, class DestAccessor>
    01695 inline void
    01696 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
    01697                                DestTraverser dul, DestAccessor dest, fftw_real norm)
    01698 {
    01699     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01700                                       norm, FFTW_REDFT00, FFTW_RODFT00);
    01701 }
    01702
    01703 template <class DestTraverser, class DestAccessor>
    01704 inline void
    01705 fourierTransformRealEO(
    01706          FFTWRealImage::const_traverser sul,
    01707          FFTWRealImage::const_traverser slr,
    01708          FFTWRealImage::Accessor src,
    01709          DestTraverser dul, DestAccessor dest, fftw_real norm)
    01710 {
    01711     int w = slr.x - sul.x;
    01712
    01713     // test for right memory layout (fftw expects a width*height fftw_real array)
    01714     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
    01715         fourierTransformRealImpl(sul, slr, dul, dest,
    01716                                  norm, FFTW_REDFT00, FFTW_RODFT00);
    01717     else
    01718         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01719                                  norm, FFTW_REDFT00, FFTW_RODFT00);
    01720 }
    01721
    01722 /********************************************************************/
    01723
    01724 template <class SrcTraverser, class SrcAccessor,
    01725           class DestTraverser, class DestAccessor>
    01726 inline void
    01727 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
    01728                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
    01729 {
    01730     fourierTransformRealOO(src.first, src.second, src.third,
    01731                                    dest.first, dest.second, norm);
    01732 }
    01733
    01734 template <class SrcTraverser, class SrcAccessor,
    01735           class DestTraverser, class DestAccessor>
    01736 inline void
    01737 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
    01738                                DestTraverser dul, DestAccessor dest, fftw_real norm)
    01739 {
    01740     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01741                                       norm, FFTW_RODFT00, FFTW_RODFT00);
    01742 }
    01743
    01744 template <class DestTraverser, class DestAccessor>
    01745 inline void
    01746 fourierTransformRealOO(
    01747          FFTWRealImage::const_traverser sul,
    01748          FFTWRealImage::const_traverser slr,
    01749          FFTWRealImage::Accessor src,
    01750          DestTraverser dul, DestAccessor dest, fftw_real norm)
    01751 {
    01752     int w = slr.x - sul.x;
    01753
    01754     // test for right memory layout (fftw expects a width*height fftw_real array)
    01755     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
    01756         fourierTransformRealImpl(sul, slr, dul, dest,
    01757                                  norm, FFTW_RODFT00, FFTW_RODFT00);
    01758     else
    01759         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
    01760                                  norm, FFTW_RODFT00, FFTW_RODFT00);
    01761 }
    01762
    01763 /*******************************************************************/
    01764
    01765 template <class SrcTraverser, class SrcAccessor,
    01766           class DestTraverser, class DestAccessor>
    01767 void
    01768 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
    01769                                   DestTraverser dul, DestAccessor dest,
    01770                                   fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
    01771 {
    01772     FFTWRealImage workImage(slr - sul);
    01773     copyImage(srcIterRange(sul, slr, src), destImage(workImage));
    01774     FFTWRealImage const & cworkImage = workImage;
    01775     fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
    01776                              dul, dest, norm, kindx, kindy);
    01777 }
    01778
    01779
    01780 template <class DestTraverser, class DestAccessor>
    01781 void
    01782 fourierTransformRealImpl(
    01783          FFTWRealImage::const_traverser sul,
    01784          FFTWRealImage::const_traverser slr,
    01785          DestTraverser dul, DestAccessor dest,
    01786          fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
    01787 {
    01788     int w = slr.x - sul.x;
    01789     int h = slr.y - sul.y;
    01790     BasicImage<fftw_real> res(w, h);
    01791
    01792     fftw_plan plan = fftw_plan_r2r_2d(h, w,
    01793                          (fftw_real *)&(*sul), (fftw_real *)res.begin(),
    01794                          kindy, kindx, FFTW_ESTIMATE);
    01795     fftw_execute(plan);
    01796     fftw_destroy_plan(plan);
    01797
    01798     if(norm != 1.0)
    01799         transformImage(srcImageRange(res), destIter(dul, dest),
    01800                        std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
    01801     else
    01802         copyImage(srcImageRange(res), destIter(dul, dest));
    01803 }
    01804
    01805
    01806 //@}
    01807
    01808 } // namespace vigra
    01809
    01810 #endif // VIGRA_FFTW3_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)