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

vigra/fftw3.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.7.1 (3 Dec 2010)