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

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