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

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

html generated using doxygen and Python
VIGRA 1.3.3 (18 Aug 2005)