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