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

details vigra/fftw3.hxx VIGRA

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.5.0, Dec 07 2006 )                                    */
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 /*        koethe@informatik.uni-hamburg.de          or                  */
00012 /*        vigra@kogs1.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 3. There is also an
00870     \link FourierTransformFFTW2 interface to the older version FFTW 2\endlink
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 template <class SrcImageIterator, class SrcAccessor,
01012           class DestImageIterator, class DestAccessor>
01013 void moveDCToCenter(SrcImageIterator src_upperleft,
01014                                SrcImageIterator src_lowerright, SrcAccessor sa,
01015                                DestImageIterator dest_upperleft, DestAccessor da)
01016 {
01017     int w= src_lowerright.x - src_upperleft.x;
01018     int h= src_lowerright.y - src_upperleft.y;
01019     int w1 = w/2;
01020     int h1 = h/2;
01021     int w2 = (w+1)/2;
01022     int h2 = (h+1)/2;
01023 
01024     // 2. Quadrant  zum 4.
01025     copyImage(srcIterRange(src_upperleft,
01026                            src_upperleft  + Diff2D(w2, h2), sa),
01027               destIter    (dest_upperleft + Diff2D(w1, h1), da));
01028 
01029     // 4. Quadrant zum 2.
01030     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
01031                            src_lowerright, sa),
01032               destIter    (dest_upperleft, da));
01033 
01034     // 1. Quadrant zum 3.
01035     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
01036                            src_upperleft  + Diff2D(w,  h2), sa),
01037               destIter    (dest_upperleft + Diff2D(0,  h1), da));
01038 
01039     // 3. Quadrant zum 1.
01040     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
01041                            src_upperleft  + Diff2D(w2, h), sa),
01042               destIter    (dest_upperleft + Diff2D(w1, 0), da));
01043 }
01044 
01045 template <class SrcImageIterator, class SrcAccessor,
01046           class DestImageIterator, class DestAccessor>
01047 inline void moveDCToCenter(
01048     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01049     pair<DestImageIterator, DestAccessor> dest)
01050 {
01051     moveDCToCenter(src.first, src.second, src.third,
01052                           dest.first, dest.second);
01053 }
01054 
01055 /********************************************************/
01056 /*                                                      */
01057 /*                   moveDCToUpperLeft                  */
01058 /*                                                      */
01059 /********************************************************/
01060 
01061 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
01062           in the image's upper left.
01063 
01064      This function is the inversion of \ref moveDCToCenter(). See there
01065      for declarations and a usage example.
01066 
01067     <b> Declarations:</b>
01068 
01069     pass arguments explicitly:
01070     \code
01071     namespace vigra {
01072         template <class SrcImageIterator, class SrcAccessor,
01073           class DestImageIterator, class DestAccessor>
01074         void moveDCToUpperLeft(SrcImageIterator src_upperleft,
01075                                SrcImageIterator src_lowerright, SrcAccessor sa,
01076                                DestImageIterator dest_upperleft, DestAccessor da);
01077     }
01078     \endcode
01079 
01080 
01081     use argument objects in conjunction with \ref ArgumentObjectFactories:
01082     \code
01083     namespace vigra {
01084         template <class SrcImageIterator, class SrcAccessor,
01085                   class DestImageIterator, class DestAccessor>
01086         void moveDCToUpperLeft(
01087             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01088             pair<DestImageIterator, DestAccessor> dest);
01089     }
01090     \endcode
01091 */
01092 template <class SrcImageIterator, class SrcAccessor,
01093           class DestImageIterator, class DestAccessor>
01094 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
01095                                SrcImageIterator src_lowerright, SrcAccessor sa,
01096                                DestImageIterator dest_upperleft, DestAccessor da)
01097 {
01098     int w= src_lowerright.x - src_upperleft.x;
01099     int h= src_lowerright.y - src_upperleft.y;
01100     int w2 = w/2;
01101     int h2 = h/2;
01102     int w1 = (w+1)/2;
01103     int h1 = (h+1)/2;
01104 
01105     // 2. Quadrant  zum 4.
01106     copyImage(srcIterRange(src_upperleft,
01107                            src_upperleft  + Diff2D(w2, h2), sa),
01108               destIter    (dest_upperleft + Diff2D(w1, h1), da));
01109 
01110     // 4. Quadrant zum 2.
01111     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
01112                            src_lowerright, sa),
01113               destIter    (dest_upperleft, da));
01114 
01115     // 1. Quadrant zum 3.
01116     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
01117                            src_upperleft  + Diff2D(w,  h2), sa),
01118               destIter    (dest_upperleft + Diff2D(0,  h1), da));
01119 
01120     // 3. Quadrant zum 1.
01121     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
01122                            src_upperleft  + Diff2D(w2, h), sa),
01123               destIter    (dest_upperleft + Diff2D(w1, 0), da));
01124 }
01125 
01126 template <class SrcImageIterator, class SrcAccessor,
01127           class DestImageIterator, class DestAccessor>
01128 inline void moveDCToUpperLeft(
01129     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01130     pair<DestImageIterator, DestAccessor> dest)
01131 {
01132     moveDCToUpperLeft(src.first, src.second, src.third,
01133                                           dest.first, dest.second);
01134 }
01135 
01136 namespace detail {
01137 
01138 template <class T>
01139 void 
01140 fourierTransformImpl(FFTWComplexImage::const_traverser sul,
01141                      FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01142                      FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest, T sign)
01143 {
01144     int w = slr.x - sul.x;
01145     int h = slr.y - sul.y;
01146 
01147     FFTWComplexImage sworkImage, dworkImage;
01148     
01149     fftw_complex * srcPtr = (fftw_complex *)(&*sul);
01150     fftw_complex * destPtr = (fftw_complex *)(&*dul);
01151     
01152     // test for right memory layout (fftw expects a 2*width*height floats array)
01153     if (&(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1))))
01154     {
01155         sworkImage.resize(w, h);
01156         copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
01157         srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft()));
01158     }
01159     if (&(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
01160     {
01161         dworkImage.resize(w, h);
01162         destPtr = (fftw_complex *)(&(*dworkImage.upperLeft()));
01163     }
01164 
01165     fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
01166     fftw_execute(plan);
01167     fftw_destroy_plan(plan);
01168     
01169     if (&(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
01170     {
01171         copyImage(srcImageRange(dworkImage), destIter(dul, dest));
01172     }
01173 }
01174 
01175 } // namespace detail
01176 
01177 /********************************************************/
01178 /*                                                      */
01179 /*                    fourierTrasform                   */
01180 /*                                                      */
01181 /********************************************************/
01182 
01183 /** \brief Compute forward and inverse Fourier transforms.
01184 
01185     In the forward direction, the input image may be scalar or complex, and the output image
01186     is always complex. In the inverse direction, both input and output must be complex.    
01187     
01188     <b> Declarations:</b>
01189 
01190     pass arguments explicitly:
01191     \code
01192     namespace vigra {
01193         template <class SrcImageIterator, class SrcAccessor>
01194         void fourierTransform(SrcImageIterator srcUpperLeft,
01195                               SrcImageIterator srcLowerRight, SrcAccessor src,
01196                               FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest);
01197 
01198         void 
01199         fourierTransformInverse(FFTWComplexImage::const_traverser sul,
01200                                 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01201                                 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
01202     }
01203     \endcode
01204 
01205     use argument objects in conjunction with \ref ArgumentObjectFactories:
01206     \code
01207     namespace vigra {
01208         template <class SrcImageIterator, class SrcAccessor>
01209         void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01210                               pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
01211 
01212         void 
01213         fourierTransformInverse(triple<FFTWComplexImage::const_traverser, 
01214                                        FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
01215                                 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
01216     }
01217     \endcode
01218 
01219     <b> Usage:</b>
01220 
01221     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
01222     Namespace: vigra
01223 
01224     \code
01225     // compute complex Fourier transform of a real image
01226     vigra::DImage src(w, h);
01227     vigra::FFTWComplexImage fourier(w, h);
01228     
01229     fourierTransform(srcImageRange(src), destImage(fourier));
01230 
01231     // compute inverse Fourier transform
01232     // note that both source and destination image must be of type vigra::FFTWComplexImage
01233     vigra::FFTWComplexImage inverseFourier(w, h);
01234 
01235     fourierTransform(srcImageRange(fourier), destImage(inverseFourier));
01236     \endcode
01237 */
01238 inline void 
01239 fourierTransform(FFTWComplexImage::const_traverser sul,
01240                  FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01241                  FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
01242 {
01243     detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
01244 }
01245 
01246 template <class SrcImageIterator, class SrcAccessor>
01247 void fourierTransform(SrcImageIterator srcUpperLeft,
01248                       SrcImageIterator srcLowerRight, SrcAccessor sa,
01249                       FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor da)
01250 {
01251     // copy real input images into a complex one...
01252     int w= srcLowerRight.x - srcUpperLeft.x;
01253     int h= srcLowerRight.y - srcUpperLeft.y;
01254 
01255     FFTWComplexImage workImage(w, h);
01256     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01257               destImage(workImage, FFTWWriteRealAccessor()));
01258 
01259     // ...and call the complex -> complex version of the algorithm
01260     FFTWComplexImage const & cworkImage = workImage;
01261     fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01262                      destUpperLeft, da);
01263 }
01264 
01265 template <class SrcImageIterator, class SrcAccessor>
01266 inline
01267 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01268                       pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
01269 {
01270     fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
01271 }
01272 
01273 inline void 
01274 fourierTransformInverse(FFTWComplexImage::const_traverser sul,
01275                         FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
01276                         FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
01277 {
01278     detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
01279 }
01280 
01281 inline void 
01282 fourierTransformInverse(triple<FFTWComplexImage::const_traverser, 
01283                                FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
01284                         pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
01285 {
01286     fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second);
01287 }
01288 
01289 /********************************************************/
01290 /*                                                      */
01291 /*                   applyFourierFilter                 */
01292 /*                                                      */
01293 /********************************************************/
01294 
01295 /** \brief Apply a filter (defined in the frequency domain) to an image.
01296 
01297     After transferring the image into the frequency domain, it is
01298     multiplied pixel-wise with the filter and transformed back. The
01299     result is put into the given destination image which must have the right size. 
01300     The result will be normalized to compensate for the two FFTs. 
01301     
01302     If the destination image is scalar, only the real part of the result image is 
01303     retained. In this case, you are responsible for choosing a filter image
01304     which ensures a zero imaginary part of the result (e.g. use a real, even symmetric 
01305     filter image, or a purely imaginary, odd symmetric on).
01306 
01307     The DC entry of the filter must be in the upper left, which is the
01308     position where FFTW expects it (see \ref moveDCToUpperLeft()).
01309 
01310     <b> Declarations:</b>
01311 
01312     pass arguments explicitly:
01313     \code
01314     namespace vigra {
01315         template <class SrcImageIterator, class SrcAccessor,
01316                   class FilterImageIterator, class FilterAccessor,
01317                   class DestImageIterator, class DestAccessor>
01318         void applyFourierFilter(SrcImageIterator srcUpperLeft,
01319                                 SrcImageIterator srcLowerRight, SrcAccessor sa,
01320                                 FilterImageIterator filterUpperLeft, FilterAccessor fa,
01321                                 DestImageIterator destUpperLeft, DestAccessor da);
01322     }
01323     \endcode
01324 
01325     use argument objects in conjunction with \ref ArgumentObjectFactories:
01326     \code
01327     namespace vigra {
01328         template <class SrcImageIterator, class SrcAccessor,
01329                 class FilterImageIterator, class FilterAccessor,
01330                 class DestImageIterator, class DestAccessor>
01331         void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01332                                 pair<FilterImageIterator, FilterAccessor> filter,
01333                                 pair<DestImageIterator, DestAccessor> dest);
01334     }
01335     \endcode
01336 
01337     <b> Usage:</b>
01338 
01339     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
01340     Namespace: vigra
01341 
01342     \code
01343     // create a Gaussian filter in Fourier space
01344     vigra::FImage gaussFilter(w, h), filter(w, h);
01345     for(int y=0; y<h; ++y)
01346         for(int x=0; x<w; ++x)
01347         {
01348             xx = float(x - w / 2) / w;
01349             yy = float(y - h / 2) / h;
01350 
01351             gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
01352         }
01353 
01354     // applyFourierFilter() expects the filter's DC in the upper left
01355     moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
01356 
01357     vigra::FFTWComplexImage result(w, h);
01358 
01359     vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
01360     \endcode
01361 
01362     For inspection of the result, \ref FFTWMagnitudeAccessor might be
01363     useful. If you want to apply the same filter repeatedly, it may be more
01364     efficient to use the FFTW functions directly with FFTW plans optimized
01365     for good performance.
01366 */
01367 template <class SrcImageIterator, class SrcAccessor,
01368           class FilterImageIterator, class FilterAccessor,
01369           class DestImageIterator, class DestAccessor>
01370 void applyFourierFilter(SrcImageIterator srcUpperLeft,
01371                         SrcImageIterator srcLowerRight, SrcAccessor sa,
01372                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
01373                         DestImageIterator destUpperLeft, DestAccessor da)
01374 {
01375     // copy real input images into a complex one...
01376     int w= srcLowerRight.x - srcUpperLeft.x;
01377     int h= srcLowerRight.y - srcUpperLeft.y;
01378 
01379     FFTWComplexImage workImage(w, h);
01380     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01381               destImage(workImage, FFTWWriteRealAccessor()));
01382 
01383     // ...and call the impl
01384     FFTWComplexImage const & cworkImage = workImage;
01385     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01386                            filterUpperLeft, fa,
01387                            destUpperLeft, da);
01388 }
01389 
01390 template <class FilterImageIterator, class FilterAccessor,
01391           class DestImageIterator, class DestAccessor>
01392 inline
01393 void applyFourierFilter(
01394     FFTWComplexImage::const_traverser srcUpperLeft,
01395     FFTWComplexImage::const_traverser srcLowerRight,
01396     FFTWComplexImage::ConstAccessor sa,
01397     FilterImageIterator filterUpperLeft, FilterAccessor fa,
01398     DestImageIterator destUpperLeft, DestAccessor da)
01399 {
01400     int w = srcLowerRight.x - srcUpperLeft.x;
01401     int h = srcLowerRight.y - srcUpperLeft.y;
01402 
01403     // test for right memory layout (fftw expects a 2*width*height floats array)
01404     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
01405         applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
01406                                filterUpperLeft, fa,
01407                                destUpperLeft, da);
01408     else
01409     {
01410         FFTWComplexImage workImage(w, h);
01411         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01412                   destImage(workImage));
01413 
01414         FFTWComplexImage const & cworkImage = workImage;
01415         applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01416                                filterUpperLeft, fa,
01417                                destUpperLeft, da);
01418     }
01419 }
01420 
01421 template <class SrcImageIterator, class SrcAccessor,
01422           class FilterImageIterator, class FilterAccessor,
01423           class DestImageIterator, class DestAccessor>
01424 inline
01425 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01426                         pair<FilterImageIterator, FilterAccessor> filter,
01427                         pair<DestImageIterator, DestAccessor> dest)
01428 {
01429     applyFourierFilter(src.first, src.second, src.third,
01430                        filter.first, filter.second,
01431                        dest.first, dest.second);
01432 }
01433 
01434 template <class FilterImageIterator, class FilterAccessor,
01435           class DestImageIterator, class DestAccessor>
01436 void applyFourierFilterImpl(
01437     FFTWComplexImage::const_traverser srcUpperLeft,
01438     FFTWComplexImage::const_traverser srcLowerRight,
01439     FFTWComplexImage::ConstAccessor sa,
01440     FilterImageIterator filterUpperLeft, FilterAccessor fa,
01441     DestImageIterator destUpperLeft, DestAccessor da)
01442 {
01443     int w = srcLowerRight.x - srcUpperLeft.x;
01444     int h = srcLowerRight.y - srcUpperLeft.y;
01445 
01446     FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
01447 
01448     // FFT from srcImage to complexResultImg
01449     fftw_plan forwardPlan=
01450         fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
01451                                (fftw_complex *)complexResultImg.begin(),
01452                                FFTW_FORWARD, FFTW_ESTIMATE );
01453     fftw_execute(forwardPlan);
01454     fftw_destroy_plan(forwardPlan);
01455 
01456     // convolve in freq. domain (in complexResultImg)
01457     combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
01458                      destImage(complexResultImg), std::multiplies<FFTWComplex>());
01459 
01460     // FFT back into spatial domain (inplace in complexResultImg)
01461     fftw_plan backwardPlan=
01462         fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
01463                                (fftw_complex *)complexResultImg.begin(),
01464                                FFTW_BACKWARD, FFTW_ESTIMATE);
01465     fftw_execute(backwardPlan);
01466     fftw_destroy_plan(backwardPlan);
01467 
01468     typedef typename
01469         NumericTraits<typename DestAccessor::value_type>::isScalar
01470         isScalarResult;
01471 
01472     // normalization (after FFTs), maybe stripping imaginary part
01473     applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
01474                                         isScalarResult());
01475 }
01476 
01477 template <class DestImageIterator, class DestAccessor>
01478 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
01479                                          DestImageIterator destUpperLeft,
01480                                          DestAccessor da,
01481                                          VigraFalseType)
01482 {
01483     double normFactor= 1.0/(srcImage.width() * srcImage.height());
01484 
01485     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
01486     {
01487         DestImageIterator dIt= destUpperLeft;
01488         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
01489         {
01490             da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
01491             da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
01492         }
01493     }
01494 }
01495 
01496 inline
01497 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
01498         FFTWComplexImage::traverser destUpperLeft,
01499         FFTWComplexImage::Accessor da,
01500         VigraFalseType)
01501 {
01502     transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
01503                    linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height())));
01504 }
01505 
01506 template <class DestImageIterator, class DestAccessor>
01507 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
01508                                          DestImageIterator destUpperLeft,
01509                                          DestAccessor da,
01510                                          VigraTrueType)
01511 {
01512     double normFactor= 1.0/(srcImage.width() * srcImage.height());
01513 
01514     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
01515     {
01516         DestImageIterator dIt= destUpperLeft;
01517         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
01518             da.set(srcImage(x, y).re()*normFactor, dIt);
01519     }
01520 }
01521 
01522 /**********************************************************/
01523 /*                                                        */
01524 /*                applyFourierFilterFamily                */
01525 /*                                                        */
01526 /**********************************************************/
01527 
01528 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
01529 
01530     This provides the same functionality as \ref applyFourierFilter(),
01531     but applying several filters at once allows to avoid
01532     repeated Fourier transforms of the source image.
01533 
01534     Filters and result images must be stored in \ref vigra::ImageArray data
01535     structures. In contrast to \ref applyFourierFilter(), this function adjusts
01536     the size of the result images and the the length of the array.
01537 
01538     <b> Declarations:</b>
01539 
01540     pass arguments explicitly:
01541     \code
01542     namespace vigra {
01543         template <class SrcImageIterator, class SrcAccessor, class FilterType>
01544         void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
01545                                       SrcImageIterator srcLowerRight, SrcAccessor sa,
01546                                       const ImageArray<FilterType> &filters,
01547                                       ImageArray<FFTWComplexImage> &results)
01548     }
01549     \endcode
01550 
01551     use argument objects in conjunction with \ref ArgumentObjectFactories:
01552     \code
01553     namespace vigra {
01554         template <class SrcImageIterator, class SrcAccessor, class FilterType>
01555         inline
01556         void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01557                                       const ImageArray<FilterType> &filters,
01558                                       ImageArray<FFTWComplexImage> &results)
01559     }
01560     \endcode
01561 
01562     <b> Usage:</b>
01563 
01564     <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
01565     Namespace: vigra
01566 
01567     \code
01568     // assuming the presence of a real-valued image named "image" to
01569     // be filtered in this example
01570 
01571     vigra::ImageArray<vigra::FImage> filters(16, image.size());
01572 
01573     for (int i=0; i<filters.size(); i++)
01574          // create some meaningful filters here
01575          createMyFilterOfScale(i, destImage(filters[i]));
01576 
01577     vigra::ImageArray<vigra::FFTWComplexImage> results();
01578 
01579     vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
01580     \endcode
01581 */
01582 template <class SrcImageIterator, class SrcAccessor,
01583           class FilterType, class DestImage>
01584 inline
01585 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01586                               const ImageArray<FilterType> &filters,
01587                               ImageArray<DestImage> &results)
01588 {
01589     applyFourierFilterFamily(src.first, src.second, src.third,
01590                              filters, results);
01591 }
01592 
01593 template <class SrcImageIterator, class SrcAccessor,
01594           class FilterType, class DestImage>
01595 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
01596                               SrcImageIterator srcLowerRight, SrcAccessor sa,
01597                               const ImageArray<FilterType> &filters,
01598                               ImageArray<DestImage> &results)
01599 {
01600     int w= srcLowerRight.x - srcUpperLeft.x;
01601     int h= srcLowerRight.y - srcUpperLeft.y;
01602 
01603     FFTWComplexImage workImage(w, h);
01604     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01605               destImage(workImage, FFTWWriteRealAccessor()));
01606 
01607     FFTWComplexImage const & cworkImage = workImage;
01608     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01609                                  filters, results);
01610 }
01611 
01612 template <class FilterType, class DestImage>
01613 inline
01614 void applyFourierFilterFamily(
01615     FFTWComplexImage::const_traverser srcUpperLeft,
01616     FFTWComplexImage::const_traverser srcLowerRight,
01617     FFTWComplexImage::ConstAccessor sa,
01618     const ImageArray<FilterType> &filters,
01619     ImageArray<DestImage> &results)
01620 {
01621     int w= srcLowerRight.x - srcUpperLeft.x;
01622 
01623     // test for right memory layout (fftw expects a 2*width*height floats array)
01624     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
01625         applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
01626                                      filters, results);
01627     else
01628     {
01629         int h = srcLowerRight.y - srcUpperLeft.y;
01630         FFTWComplexImage workImage(w, h);
01631         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01632                   destImage(workImage));
01633 
01634         FFTWComplexImage const & cworkImage = workImage;
01635         applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01636                                      filters, results);
01637     }
01638 }
01639 
01640 template <class FilterType, class DestImage>
01641 void applyFourierFilterFamilyImpl(
01642     FFTWComplexImage::const_traverser srcUpperLeft,
01643     FFTWComplexImage::const_traverser srcLowerRight,
01644     FFTWComplexImage::ConstAccessor sa,
01645     const ImageArray<FilterType> &filters,
01646     ImageArray<DestImage> &results)
01647 {
01648     // make sure the filter images have the right dimensions
01649     vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
01650                        "applyFourierFilterFamily called with src image size != filters.imageSize()!");
01651 
01652     // make sure the result image array has the right dimensions
01653     results.resize(filters.size());
01654     results.resizeImages(filters.imageSize());
01655 
01656     // FFT from srcImage to freqImage
01657     int w= srcLowerRight.x - srcUpperLeft.x;
01658     int h= srcLowerRight.y - srcUpperLeft.y;
01659 
01660     FFTWComplexImage freqImage(w, h);
01661     FFTWComplexImage result(w, h);
01662 
01663     fftw_plan forwardPlan=
01664         fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
01665                                (fftw_complex *)freqImage.begin(),
01666                                FFTW_FORWARD, FFTW_ESTIMATE );
01667     fftw_execute(forwardPlan);
01668     fftw_destroy_plan(forwardPlan);
01669 
01670     fftw_plan backwardPlan=
01671         fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
01672                                (fftw_complex *)result.begin(),
01673                                FFTW_BACKWARD, FFTW_ESTIMATE );
01674     typedef typename
01675         NumericTraits<typename DestImage::Accessor::value_type>::isScalar
01676         isScalarResult;
01677 
01678     // convolve with filters in freq. domain
01679     for (unsigned int i= 0;  i < filters.size(); i++)
01680     {
01681         combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
01682                          destImage(result), std::multiplies<FFTWComplex>());
01683 
01684         // FFT back into spatial domain (inplace in destImage)
01685         fftw_execute(backwardPlan);
01686 
01687         // normalization (after FFTs), maybe stripping imaginary part
01688         applyFourierFilterImplNormalization(result,
01689                                             results[i].upperLeft(), results[i].accessor(),
01690                                             isScalarResult());
01691     }
01692     fftw_destroy_plan(backwardPlan);
01693 }
01694 
01695 /********************************************************/
01696 /*                                                      */
01697 /*                fourierTransformReal                  */
01698 /*                                                      */
01699 /********************************************************/
01700 
01701 /** \brief Real Fourier transforms for even and odd boundary conditions
01702            (aka. cosine and sine transforms).
01703 
01704     
01705     If the image is real and has even symmetry, its Fourier transform
01706     is also real and has even symmetry. The Fourier transform of a real image with odd
01707     symmetry is imaginary and has odd symmetry. In either case, only about a quarter
01708     of the pixels need to be stored because the rest can be calculated from the symmetry
01709     properties. This is especially useful, if the original image is implicitly assumed
01710     to have reflective or anti-reflective boundary conditions. Then the "negative"
01711     pixel locations are defined as
01712 
01713     \code
01714     even (reflective boundary conditions):      f[-x] = f[x]     (x = 1,...,N-1)
01715     odd (anti-reflective boundary conditions):  f[-1] = 0
01716                                                 f[-x] = -f[x-2]  (x = 2,...,N-1)
01717     \endcode
01718     
01719     end similar at the other boundary (see the FFTW documentation for details). 
01720     This has the advantage that more efficient Fourier transforms that use only
01721     real numbers can be implemented. These are also known as cosine and sine transforms
01722     respectively. 
01723     
01724     If you use the odd transform it is important to note that in the Fourier domain,
01725     the DC component is always zero and is therefore dropped from the data structure.
01726     This means that index 0 in an odd symmetric Fourier domain image refers to
01727     the <i>first</i> harmonic. This is especially important if an image is first 
01728     cosine transformed (even symmetry), then in the Fourier domain multiplied 
01729     with an odd symmetric filter (e.g. a first derivative) and finally transformed
01730     back to the spatial domain with a sine transform (odd symmetric). For this to work
01731     properly the image must be shifted left or up by one pixel (depending on whether
01732     the x- or y-axis is odd symmetric) before the inverse transform can be applied.
01733     (see example below).
01734     
01735     The real Fourier transform functions are named <tt>fourierTransformReal??</tt>
01736     where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating
01737     whether the x- and y-axis is to be transformed using even or odd symmetry.
01738     The same functions can be used for both the forward and inverse transforms,
01739     only the normalization changes. For signal processing, the following 
01740     normalization factors are most appropriate:
01741     
01742     \code
01743                           forward             inverse
01744     ------------------------------------------------------------
01745     X even, Y even           1.0         4.0 * (w-1) * (h-1)
01746     X even, Y odd           -1.0        -4.0 * (w-1) * (h+1)
01747     X odd,  Y even          -1.0        -4.0 * (w+1) * (h-1)
01748     X odd,  Y odd            1.0         4.0 * (w+1) * (h+1)
01749     \endcode
01750 
01751     where <tt>w</tt> and <tt>h</tt> denote the image width and height.
01752 
01753     <b> Declarations:</b>
01754 
01755     pass arguments explicitly:
01756     \code
01757     namespace vigra {
01758         template <class SrcTraverser, class SrcAccessor,
01759                   class DestTraverser, class DestAccessor>
01760         void
01761         fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
01762                                DestTraverser dul, DestAccessor dest, fftw_real norm);
01763                                
01764         fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
01765     }
01766     \endcode
01767 
01768 
01769     use argument objects in conjunction with \ref ArgumentObjectFactories:
01770     \code
01771     namespace vigra {
01772         template <class SrcTraverser, class SrcAccessor,
01773                   class DestTraverser, class DestAccessor>
01774         void
01775         fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
01776                                pair<DestTraverser, DestAccessor> dest, fftw_real norm);
01777                                
01778         fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
01779     }
01780     \endcode
01781 
01782     <b> Usage:</b>
01783 
01784         <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br>
01785         Namespace: vigra
01786 
01787     \code
01788     vigra::FImage spatial(width,height), fourier(width,height);
01789     ... // fill image with data
01790 
01791     // forward cosine transform == reflective boundary conditions
01792     fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
01793 
01794     // multiply with a first derivative of Gaussian in x-direction
01795     for(int y = 0; y < height; ++y)
01796     {
01797         for(int x = 1; x < width; ++x)
01798         {
01799             double dx = x * M_PI / (width - 1);
01800             double dy = y * M_PI / (height - 1);
01801             fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
01802         }
01803         fourier(width-1, y) = 0.0;
01804     }
01805     
01806     // inverse transform -- odd symmetry in x-direction, even in y,
01807     //                      due to symmetry of the filter
01808     fourierTransformRealOE(srcImageRange(fourier), destImage(spatial), 
01809                            (fftw_real)-4.0 * (width+1) * (height-1));
01810     \endcode
01811 */
01812 template <class SrcTraverser, class SrcAccessor,
01813           class DestTraverser, class DestAccessor>
01814 inline void
01815 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
01816                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
01817 {
01818     fourierTransformRealEE(src.first, src.second, src.third,
01819                                    dest.first, dest.second, norm);
01820 }
01821 
01822 template <class SrcTraverser, class SrcAccessor,
01823           class DestTraverser, class DestAccessor>
01824 inline void
01825 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
01826                                DestTraverser dul, DestAccessor dest, fftw_real norm)
01827 {
01828     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01829                                       norm, FFTW_REDFT00, FFTW_REDFT00);
01830 }
01831 
01832 template <class DestTraverser, class DestAccessor>
01833 inline void
01834 fourierTransformRealEE(
01835          FFTWRealImage::const_traverser sul,
01836          FFTWRealImage::const_traverser slr,
01837          FFTWRealImage::Accessor src,
01838          DestTraverser dul, DestAccessor dest, fftw_real norm)
01839 {
01840     int w = slr.x - sul.x;
01841 
01842     // test for right memory layout (fftw expects a width*height fftw_real array)
01843     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
01844         fourierTransformRealImpl(sul, slr, dul, dest,
01845                                  norm, FFTW_REDFT00, FFTW_REDFT00);
01846     else
01847         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01848                                  norm, FFTW_REDFT00, FFTW_REDFT00);
01849 }
01850 
01851 /********************************************************************/
01852 
01853 template <class SrcTraverser, class SrcAccessor,
01854           class DestTraverser, class DestAccessor>
01855 inline void
01856 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
01857                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
01858 {
01859     fourierTransformRealOE(src.first, src.second, src.third,
01860                                    dest.first, dest.second, norm);
01861 }
01862 
01863 template <class SrcTraverser, class SrcAccessor,
01864           class DestTraverser, class DestAccessor>
01865 inline void
01866 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
01867                                DestTraverser dul, DestAccessor dest, fftw_real norm)
01868 {
01869     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01870                                       norm, FFTW_RODFT00, FFTW_REDFT00);
01871 }
01872 
01873 template <class DestTraverser, class DestAccessor>
01874 inline void
01875 fourierTransformRealOE(
01876          FFTWRealImage::const_traverser sul,
01877          FFTWRealImage::const_traverser slr,
01878          FFTWRealImage::Accessor src,
01879          DestTraverser dul, DestAccessor dest, fftw_real norm)
01880 {
01881     int w = slr.x - sul.x;
01882 
01883     // test for right memory layout (fftw expects a width*height fftw_real array)
01884     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
01885         fourierTransformRealImpl(sul, slr, dul, dest,
01886                                  norm, FFTW_RODFT00, FFTW_REDFT00);
01887     else
01888         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01889                                  norm, FFTW_RODFT00, FFTW_REDFT00);
01890 }
01891 
01892 /********************************************************************/
01893 
01894 template <class SrcTraverser, class SrcAccessor,
01895           class DestTraverser, class DestAccessor>
01896 inline void
01897 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
01898                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
01899 {
01900     fourierTransformRealEO(src.first, src.second, src.third,
01901                                    dest.first, dest.second, norm);
01902 }
01903 
01904 template <class SrcTraverser, class SrcAccessor,
01905           class DestTraverser, class DestAccessor>
01906 inline void
01907 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
01908                                DestTraverser dul, DestAccessor dest, fftw_real norm)
01909 {
01910     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01911                                       norm, FFTW_REDFT00, FFTW_RODFT00);
01912 }
01913 
01914 template <class DestTraverser, class DestAccessor>
01915 inline void
01916 fourierTransformRealEO(
01917          FFTWRealImage::const_traverser sul,
01918          FFTWRealImage::const_traverser slr,
01919          FFTWRealImage::Accessor src,
01920          DestTraverser dul, DestAccessor dest, fftw_real norm)
01921 {
01922     int w = slr.x - sul.x;
01923 
01924     // test for right memory layout (fftw expects a width*height fftw_real array)
01925     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
01926         fourierTransformRealImpl(sul, slr, dul, dest,
01927                                  norm, FFTW_REDFT00, FFTW_RODFT00);
01928     else
01929         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01930                                  norm, FFTW_REDFT00, FFTW_RODFT00);
01931 }
01932 
01933 /********************************************************************/
01934 
01935 template <class SrcTraverser, class SrcAccessor,
01936           class DestTraverser, class DestAccessor>
01937 inline void
01938 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
01939                                pair<DestTraverser, DestAccessor> dest, fftw_real norm)
01940 {
01941     fourierTransformRealOO(src.first, src.second, src.third,
01942                                    dest.first, dest.second, norm);
01943 }
01944 
01945 template <class SrcTraverser, class SrcAccessor,
01946           class DestTraverser, class DestAccessor>
01947 inline void
01948 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
01949                                DestTraverser dul, DestAccessor dest, fftw_real norm)
01950 {
01951     fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01952                                       norm, FFTW_RODFT00, FFTW_RODFT00);
01953 }
01954 
01955 template <class DestTraverser, class DestAccessor>
01956 inline void
01957 fourierTransformRealOO(
01958          FFTWRealImage::const_traverser sul,
01959          FFTWRealImage::const_traverser slr,
01960          FFTWRealImage::Accessor src,
01961          DestTraverser dul, DestAccessor dest, fftw_real norm)
01962 {
01963     int w = slr.x - sul.x;
01964 
01965     // test for right memory layout (fftw expects a width*height fftw_real array)
01966     if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
01967         fourierTransformRealImpl(sul, slr, dul, dest,
01968                                  norm, FFTW_RODFT00, FFTW_RODFT00);
01969     else
01970         fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
01971                                  norm, FFTW_RODFT00, FFTW_RODFT00);
01972 }
01973 
01974 /*******************************************************************/
01975 
01976 template <class SrcTraverser, class SrcAccessor,
01977           class DestTraverser, class DestAccessor>
01978 void
01979 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
01980                                   DestTraverser dul, DestAccessor dest,
01981                                   fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
01982 {
01983     FFTWRealImage workImage(slr - sul);
01984     copyImage(srcIterRange(sul, slr, src), destImage(workImage));
01985     FFTWRealImage const & cworkImage = workImage;
01986     fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
01987                              dul, dest, norm, kindx, kindy);
01988 }
01989 
01990 
01991 template <class DestTraverser, class DestAccessor>
01992 void
01993 fourierTransformRealImpl(
01994          FFTWRealImage::const_traverser sul,
01995          FFTWRealImage::const_traverser slr,
01996          DestTraverser dul, DestAccessor dest,
01997          fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
01998 {
01999     int w = slr.x - sul.x;
02000     int h = slr.y - sul.y;
02001     BasicImage<fftw_real> res(w, h);
02002     
02003     fftw_plan plan = fftw_plan_r2r_2d(h, w, 
02004                          (fftw_real *)&(*sul), (fftw_real *)res.begin(),
02005                          kindy, kindx, FFTW_ESTIMATE);
02006     fftw_execute(plan);
02007     fftw_destroy_plan(plan);
02008 
02009     if(norm != 1.0)
02010         transformImage(srcImageRange(res), destIter(dul, dest),
02011                        std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
02012     else
02013         copyImage(srcImageRange(res), destIter(dul, dest));
02014 }
02015 
02016 
02017 //@}
02018 
02019 } // namespace vigra
02020 
02021 #endif // VIGRA_FFTW3_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)