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

vigra/fftw3.hxx

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)