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

details vigra/fftw.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.2.0, Aug 07 2003 )                                    */
00008 /*    You may use, modify, and distribute this software according       */
00009 /*    to the terms stated in the LICENSE file included in               */
00010 /*    the VIGRA distribution.                                           */
00011 /*                                                                      */
00012 /*    The VIGRA Website is                                              */
00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00014 /*    Please direct questions, bug reports, and contributions to        */
00015 /*        koethe@informatik.uni-hamburg.de                              */
00016 /*                                                                      */
00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
00020 /*                                                                      */
00021 /************************************************************************/
00022 
00023 #ifndef VIGRA_FFTW_HXX
00024 #define VIGRA_FFTW_HXX
00025 
00026 #include <cmath>
00027 #include <functional>
00028 #include "vigra/stdimage.hxx"
00029 #include "vigra/copyimage.hxx"
00030 #include "vigra/combineimages.hxx"
00031 #include "vigra/numerictraits.hxx"
00032 #include "vigra/imagecontainer.hxx"
00033 #include <fftw.h>
00034 
00035 namespace vigra {
00036 
00037 /********************************************************/
00038 /*                                                      */
00039 /*                    FFTWComplex                       */
00040 /*                                                      */
00041 /********************************************************/
00042 
00043 /** \brief Wrapper class for the FFTW type '<TT>fftw_complex</TT>'.
00044 
00045     This class provides constructors and other member functions
00046     for the C struct '<TT>fftw_complex</TT>'. This struct is the basic
00047     pixel type of the <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a>
00048     library. It inherits the data members '<TT>re</TT>' and '<TT>im</TT>'
00049     that denote the real and imaginary part of the number. In addition it
00050     defines transformations to polar coordinates,
00051     as well as \ref FFTWComplexOperators "arithmetic operators" and
00052     \ref FFTWComplexAccessors "accessors".
00053 
00054     FFTWComplex implements the concepts \ref AlgebraicField and
00055     \ref DivisionAlgebra.
00056 
00057     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00058     Namespace: vigra
00059 */
00060 class FFTWComplex
00061 : public fftw_complex
00062 {
00063   public:
00064         /** The complex' component type, as defined in '<TT>fftw.h</TT>'
00065         */
00066     typedef fftw_real value_type;
00067 
00068         /** reference type (result of operator[])
00069         */
00070     typedef fftw_real & reference;
00071 
00072         /** const reference type (result of operator[] const)
00073         */
00074     typedef fftw_real const & const_reference;
00075 
00076         /** const reference type (result of operator[] const)
00077         */
00078     typedef fftw_real * iterator;
00079 
00080         /** const reference type (result of operator[] const)
00081         */
00082     typedef fftw_real const * const_iterator;
00083 
00084         /** Construct from real and imaginary part.
00085             Default: 0.
00086         */
00087     FFTWComplex(value_type const & ire = 0.0, value_type const & iim = 0.0)
00088     {
00089         re = ire;
00090         im = iim;
00091     }
00092 
00093         /** Copy constructor.
00094         */
00095     FFTWComplex(FFTWComplex const & o)
00096     : fftw_complex(o)
00097     {}
00098 
00099         /** Construct from plain <TT>fftw_complex</TT>.
00100         */
00101     FFTWComplex(fftw_complex const & o)
00102     : fftw_complex(o)
00103     {}
00104 
00105         /** Construct from TinyVector.
00106         */
00107     template <class T>
00108     FFTWComplex(TinyVector<T, 2> const & o)
00109     {
00110         re = o[0];
00111         im = o[1];
00112     }
00113 
00114         /** Assignment.
00115         */
00116     FFTWComplex& operator=(FFTWComplex const & o)
00117     {
00118         re = o.re;
00119         im = o.im;
00120         return *this;
00121     }
00122 
00123         /** Assignment.
00124         */
00125     FFTWComplex& operator=(fftw_complex const & o)
00126     {
00127         re = o.re;
00128         im = o.im;
00129         return *this;
00130     }
00131 
00132         /** Assignment.
00133         */
00134     FFTWComplex& operator=(fftw_real const & o)
00135     {
00136         re = o;
00137         im = 0.0;
00138         return *this;
00139     }
00140 
00141         /** Assignment.
00142         */
00143     template <class T>
00144     FFTWComplex& operator=(TinyVector<T, 2> const & o)
00145     {
00146         re = o[0];
00147         im = o[1];
00148         return *this;
00149     }
00150 
00151         /** Unary negation.
00152         */
00153     FFTWComplex operator-() const
00154         { return FFTWComplex(-re, -im); }
00155 
00156         /** Squared magnitude x*conj(x)
00157         */
00158     value_type squaredMagnitude() const
00159         { return c_re(*this)*c_re(*this)+c_im(*this)*c_im(*this); }
00160 
00161         /** Magnitude (length of radius vector).
00162         */
00163     value_type magnitude() const
00164         { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
00165 
00166         /** Phase angle.
00167         */
00168     value_type phase() const
00169         { return VIGRA_CSTD::atan2(c_im(*this),c_re(*this)); }
00170 
00171         /** Access components as if number were a vector.
00172         */
00173     reference operator[](int i)
00174         { return (&re)[i]; }
00175 
00176         /** Read components as if number were a vector.
00177         */
00178     const_reference operator[](int i) const
00179         { return (&re)[i]; }
00180 
00181         /** Length of complex number (always 2).
00182         */
00183     int size() const
00184         { return 2; }
00185         
00186     iterator begin()
00187         { return &re; }
00188         
00189     iterator end()
00190         { return &re + 2; }
00191         
00192     const_iterator begin() const
00193         { return &re; }
00194         
00195     const_iterator end() const
00196         { return &re + 2; }
00197 };
00198 
00199 /********************************************************/
00200 /*                                                      */
00201 /*                    FFTWComplex Traits                */
00202 /*                                                      */
00203 /********************************************************/
00204 
00205 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex
00206     The numeric and promote traits for FFTWComplex follow
00207     the general specifications for \ref NumericPromotionTraits.
00208     They are implemented as follows:
00209 
00210     \code
00211 
00212     template <>
00213     struct NumericTraits<FFTWComplex >
00214     {
00215         typedef fftw_complex Type;
00216         typedef fftw_complex Promote;
00217         typedef fftw_complex RealPromote;
00218         typedef VigraFalseType isIntegral;
00219         typedef VigraFalseType isScalar;
00220         typedef VigraFalseType isOrdered;
00221 
00222         // etc.
00223     };
00224 
00225     template <>
00226     struct PromoteTraits<FFTWComplex, FFTWComplex>
00227     {
00228         typedef FFTWComplex Promote;
00229     };
00230 
00231     template <>
00232     struct PromoteTraits<FFTWComplex, double>
00233     {
00234         typedef FFTWComplex Promote;
00235     };
00236 
00237     template <>
00238     struct PromoteTraits<double, FFTWComplex>
00239     {
00240         typedef FFTWComplex Promote;
00241     };
00242     \endcode
00243 
00244     where the type '<TT>fftw_complex</TT>' is defined in '<TT>fftw.h</TT>'
00245 
00246     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00247     Namespace: vigra
00248 */
00249 template<>
00250 struct NumericTraits<fftw_complex>
00251 {
00252     typedef fftw_complex Type;
00253     typedef fftw_complex Promote;
00254     typedef fftw_complex RealPromote;
00255     typedef VigraFalseType isIntegral;
00256     typedef VigraFalseType isScalar;
00257     typedef VigraFalseType isOrdered;
00258 
00259     static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
00260     static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
00261     static FFTWComplex nonZero() { return one(); }
00262 
00263     static const Promote & toPromote(const Type & v) { return v; }
00264     static const RealPromote & toRealPromote(const Type & v) { return v; }
00265     static const Type & fromPromote(const Promote & v) { return v; }
00266     static const Type & fromRealPromote(const RealPromote & v) { return v; }
00267 };
00268 
00269 template<>
00270 struct NumericTraits<FFTWComplex>
00271 : public NumericTraits<fftw_complex>
00272 {
00273     typedef FFTWComplex Type;
00274     typedef FFTWComplex Promote;
00275     typedef FFTWComplex RealPromote;
00276 };
00277 
00278 template <>
00279 struct PromoteTraits<fftw_complex, fftw_complex>
00280 {
00281     typedef fftw_complex Promote;
00282 };
00283 
00284 template <>
00285 struct PromoteTraits<fftw_complex, double>
00286 {
00287     typedef fftw_complex Promote;
00288 };
00289 
00290 template <>
00291 struct PromoteTraits<double, fftw_complex>
00292 {
00293     typedef fftw_complex Promote;
00294 };
00295 
00296 template <>
00297 struct PromoteTraits<FFTWComplex, FFTWComplex>
00298 {
00299     typedef FFTWComplex Promote;
00300 };
00301 
00302 template <>
00303 struct PromoteTraits<FFTWComplex, double>
00304 {
00305     typedef FFTWComplex Promote;
00306 };
00307 
00308 template <>
00309 struct PromoteTraits<double, FFTWComplex>
00310 {
00311     typedef FFTWComplex Promote;
00312 };
00313 
00314 
00315 /********************************************************/
00316 /*                                                      */
00317 /*                    FFTWComplex Operations            */
00318 /*                                                      */
00319 /********************************************************/
00320 
00321 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex
00322 
00323     \brief <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>
00324 
00325     These functions fulfill the requirements of an Algebraic Field.
00326     Return types are determined according to \ref FFTWComplexTraits.
00327 
00328     Namespace: vigra
00329     <p>
00330 
00331  */
00332 //@{
00333     /// equal
00334 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) {
00335     return c_re(a) == c_re(b) && c_im(a) == c_im(b);
00336 }
00337 
00338     /// not equal
00339 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) {
00340     return c_re(a) != c_re(b) || c_im(a) != c_im(b);
00341 }
00342 
00343     /// add-assignment
00344 inline FFTWComplex &operator +=(FFTWComplex &a, const FFTWComplex &b) {
00345     c_re(a) += c_re(b);
00346     c_im(a) += c_im(b);
00347     return a;
00348 }
00349 
00350     /// subtract-assignment
00351 inline FFTWComplex &operator -=(FFTWComplex &a, const FFTWComplex &b) {
00352     c_re(a) -= c_re(b);
00353     c_im(a) -= c_im(b);
00354     return a;
00355 }
00356 
00357     /// multiply-assignment
00358 inline FFTWComplex &operator *=(FFTWComplex &a, const FFTWComplex &b) {
00359     FFTWComplex::value_type t = c_re(a)*c_re(b)-c_im(a)*c_im(b);
00360     c_im(a) = c_re(a)*c_im(b)+c_im(a)*c_re(b);
00361     c_re(a) = t;
00362     return a;
00363 }
00364 
00365     /// divide-assignment
00366 inline FFTWComplex &operator /=(FFTWComplex &a, const FFTWComplex &b) {
00367     FFTWComplex::value_type sm = b.squaredMagnitude();
00368     FFTWComplex::value_type t = (c_re(a)*c_re(b)+c_im(a)*c_im(b))/sm;
00369     c_im(a) = (c_re(b)*c_im(a)-c_re(a)*c_im(b))/sm;
00370     c_re(a) = t;
00371     return a;
00372 }
00373 
00374     /// multiply-assignment with scalar double
00375 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) {
00376     c_re(a) *= b;
00377     c_im(a) *= b;
00378     return a;
00379 }
00380 
00381     /// divide-assignment with scalar double
00382 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) {
00383     c_re(a) /= b;
00384     c_im(a) /= b;
00385     return a;
00386 }
00387 
00388     /// addition
00389 inline FFTWComplex operator +(FFTWComplex a, const FFTWComplex &b) {
00390     a += b;
00391     return a;
00392 }
00393 
00394     /// subtraction
00395 inline FFTWComplex operator -(FFTWComplex a, const FFTWComplex &b) {
00396     a -= b;
00397     return a;
00398 }
00399 
00400     /// multiplication
00401 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) {
00402     a *= b;
00403     return a;
00404 }
00405 
00406     /// right multiplication with scalar double
00407 inline FFTWComplex operator *(FFTWComplex a, const double &b) {
00408     a *= b;
00409     return a;
00410 }
00411 
00412     /// left multiplication with scalar double
00413 inline FFTWComplex operator *(const double &a, FFTWComplex b) {
00414     b *= a;
00415     return b;
00416 }
00417 
00418     /// division
00419 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) {
00420     a /= b;
00421     return a;
00422 }
00423 
00424     /// right division with scalar double
00425 inline FFTWComplex operator /(FFTWComplex a, const double &b) {
00426     a /= b;
00427     return a;
00428 }
00429 
00430 using VIGRA_CSTD::abs;
00431 
00432     /// absolute value (= magnitude)
00433 inline FFTWComplex::value_type abs(const FFTWComplex &a)
00434 {
00435     return a.magnitude();
00436 }
00437 
00438     /// complex conjugate
00439 inline FFTWComplex conj(const FFTWComplex &a)
00440 {
00441     return FFTWComplex(a.re, -a.im);
00442 }
00443 
00444 //@}
00445 
00446 
00447 /** \addtogroup StandardImageTypes
00448 */
00449 //@{
00450 
00451 /********************************************************/
00452 /*                                                      */
00453 /*                      FFTWRealImage                   */
00454 /*                                                      */
00455 /********************************************************/
00456 
00457     /** Float (<tt>fftw_real</tt>) image.
00458         The type <tt>fftw_real</tt> (either <tt>float</tt> or <tt>double</tt>)
00459         is defined during compilation of fftw and imported into VIGRA
00460         from <tt>fftw.h</tt>. FFTWRealImage uses \ref vigra::BasicImageIterator 
00461         and \ref vigra::StandardAccessor and 
00462         their const counterparts to access the data.
00463         
00464         <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00465         Namespace: vigra
00466     */
00467 typedef BasicImage<fftw_real> FFTWRealImage;
00468 
00469 /********************************************************/
00470 /*                                                      */
00471 /*                     FFTWComplexImage                 */
00472 /*                                                      */
00473 /********************************************************/
00474 
00475 template<>
00476 struct IteratorTraits<
00477         BasicImageIterator<FFTWComplex, FFTWComplex **> >
00478 {
00479     typedef BasicImageIterator<FFTWComplex, FFTWComplex **>  Iterator;
00480     typedef Iterator                             iterator;
00481     typedef iterator::iterator_category          iterator_category;
00482     typedef iterator::value_type                 value_type;
00483     typedef iterator::reference                  reference;
00484     typedef iterator::index_reference            index_reference;
00485     typedef iterator::pointer                    pointer;
00486     typedef iterator::difference_type            difference_type;
00487     typedef iterator::row_iterator               row_iterator;
00488     typedef iterator::column_iterator            column_iterator;
00489     typedef VectorAccessor<FFTWComplex>          default_accessor;
00490     typedef VectorAccessor<FFTWComplex>          DefaultAccessor;
00491 };
00492 
00493 template<>
00494 struct IteratorTraits<
00495         ConstBasicImageIterator<FFTWComplex, FFTWComplex **> >
00496 {
00497     typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **>    Iterator;
00498     typedef Iterator                             iterator;
00499     typedef iterator::iterator_category          iterator_category;
00500     typedef iterator::value_type                 value_type;
00501     typedef iterator::reference                  reference;
00502     typedef iterator::index_reference            index_reference;
00503     typedef iterator::pointer                    pointer;
00504     typedef iterator::difference_type            difference_type;
00505     typedef iterator::row_iterator               row_iterator;
00506     typedef iterator::column_iterator            column_iterator;
00507     typedef VectorAccessor<FFTWComplex>          default_accessor;
00508     typedef VectorAccessor<FFTWComplex>          DefaultAccessor;
00509 };
00510 
00511     /** Complex (FFTWComplex) image.
00512         It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and
00513         their const counterparts to access the data.
00514 
00515         <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00516         Namespace: vigra
00517     */
00518 typedef BasicImage<FFTWComplex> FFTWComplexImage;
00519 
00520 //@}
00521 
00522 /********************************************************/
00523 /*                                                      */
00524 /*                  FFTWComplex-Accessors               */
00525 /*                                                      */
00526 /********************************************************/
00527 
00528 /** \addtogroup DataAccessors
00529 */
00530 //@{
00531 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex 
00532 
00533     Encapsulate access to pixels of type FFTWComplex
00534 */
00535 //@{
00536     /** Encapsulate access to the the real part of a complex number.
00537 
00538     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00539     Namespace: vigra
00540     */
00541 class FFTWRealAccessor
00542 {
00543   public:
00544 
00545         /// The accessor's value type.
00546     typedef fftw_real value_type;
00547 
00548         /// Read real part at iterator position.
00549     template <class ITERATOR>
00550     value_type operator()(ITERATOR const & i) const {
00551         return c_re(*i);
00552     }
00553 
00554         /// Read real part at offset from iterator position.
00555     template <class ITERATOR, class DIFFERENCE>
00556     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
00557         return c_re(i[d]);
00558     }
00559 
00560         /// Write real part at iterator position.
00561     template <class ITERATOR>
00562     void set(value_type const & v, ITERATOR const & i) const {
00563         c_re(*i)= v;
00564     }
00565 
00566         /// Write real part at offset from iterator position.
00567     template <class ITERATOR, class DIFFERENCE>
00568     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
00569         c_re(i[d])= v;
00570     }
00571 };
00572 
00573     /** Encapsulate access to the the imaginary part of a complex number.
00574 
00575     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00576     Namespace: vigra
00577     */
00578 class FFTWImaginaryAccessor
00579 {
00580   public:
00581         /// The accessor's value type.
00582     typedef fftw_real value_type;
00583 
00584         /// Read imaginary part at iterator position.
00585     template <class ITERATOR>
00586     value_type operator()(ITERATOR const & i) const {
00587         return c_im(*i);
00588     }
00589 
00590         /// Read imaginary part at offset from iterator position.
00591     template <class ITERATOR, class DIFFERENCE>
00592     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
00593         return c_im(i[d]);
00594     }
00595 
00596         /// Write imaginary part at iterator position.
00597     template <class ITERATOR>
00598     void set(value_type const & v, ITERATOR const & i) const {
00599         c_im(*i)= v;
00600     }
00601 
00602         /// Write imaginary part at offset from iterator position.
00603     template <class ITERATOR, class DIFFERENCE>
00604     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
00605         c_im(i[d])= v;
00606     }
00607 };
00608 
00609     /** Write a real number into a complex one. The imaginary part is set
00610         to 0.
00611 
00612     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00613     Namespace: vigra
00614     */
00615 class FFTWWriteRealAccessor: public FFTWRealAccessor
00616 {
00617   public:
00618         /// The accessor's value type.
00619     typedef fftw_real value_type;
00620 
00621         /** Write real number at iterator position. Set imaginary part
00622             to 0.
00623         */
00624     template <class ITERATOR>
00625     void set(value_type const & v, ITERATOR const & i) const {
00626         c_re(*i)= v;
00627         c_im(*i)= 0;
00628     }
00629 
00630         /** Write real number at offset from iterator position. Set imaginary part
00631             to 0.
00632         */
00633     template <class ITERATOR, class DIFFERENCE>
00634     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
00635         c_re(i[d])= v;
00636         c_im(i[d])= 0;
00637     }
00638 };
00639 
00640     /** Calculate magnitude of complex number on the fly.
00641 
00642     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00643     Namespace: vigra
00644     */
00645 class FFTWMagnitudeAccessor
00646 {
00647   public:
00648         /// The accessor's value type.
00649     typedef fftw_real value_type;
00650 
00651         /// Read magnitude at iterator position.
00652     template <class ITERATOR>
00653     value_type operator()(ITERATOR const & i) const {
00654         return (*i).magnitude();
00655     }
00656 
00657         /// Read magnitude at offset from iterator position.
00658     template <class ITERATOR, class DIFFERENCE>
00659     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
00660         return (i[d]).magnitude();
00661     }
00662 };
00663 
00664     /** Calculate phase of complex number on the fly.
00665 
00666     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00667     Namespace: vigra
00668     */
00669 class FFTWPhaseAccessor
00670 {
00671   public:
00672         /// The accessor's value type.
00673     typedef fftw_real value_type;
00674 
00675         /// Read phase at iterator position.
00676     template <class ITERATOR>
00677     value_type operator()(ITERATOR const & i) const {
00678         return (*i).phase();
00679     }
00680 
00681         /// Read phase at offset from iterator position.
00682     template <class ITERATOR, class DIFFERENCE>
00683     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
00684         return (i[d]).phase();
00685     }
00686 };
00687 
00688 //@}
00689 //@}
00690 
00691 /********************************************************/
00692 /*                                                      */
00693 /*                    Fourier Transform                 */
00694 /*                                                      */
00695 /********************************************************/
00696 
00697 /** \page FourierTransform Fast Fourier Transform
00698     VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
00699     Transform</a> package to perform Fourier transformations. VIGRA
00700     provides a wrapper for FFTW's complex number type (FFTWComplex),
00701     but FFTW's functions are used verbatim. If the image is stored as
00702     a FFTWComplexImage, a FFT is performed like this:
00703 
00704     \code
00705     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
00706     ... // fill image with data
00707 
00708     // create a plan for optimal performance
00709     fftwnd_plan forwardPlan=
00710         fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE );
00711 
00712     // calculate FFT
00713     fftwnd_one(forwardPlan, spatial.begin(), fourier.begin());
00714     \endcode
00715 
00716     Note that in the creation of a plan, the height must be given
00717     first. Note also that <TT>spatial.begin()</TT> may only be passed
00718     to <TT>fftwnd_one</TT> if the transform shall be applied to the
00719     entire image. When you want to retrict operation to an ROI, you
00720     create a copy of the ROI in an image of appropriate size.
00721 
00722     More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
00723 
00724     FFTW produces fourier images that have the DC component (the
00725     origin of the Fourier space) in the upper left corner. Often, one
00726     wants the origin in the center of the image, so that frequencies
00727     always increase towards the border of the image. This can be
00728     achieved by calling \ref moveDCToCenter(). The inverse
00729     transformation is done by \ref moveDCToUpperLeft().
00730 
00731     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00732     Namespace: vigra
00733 */
00734 
00735 /** \addtogroup FourierHelpers Helper Functions for FFT
00736     Rearrange quadrants in a Fourier image or apply filters on images
00737     using the FFTW.
00738 */
00739 //@{
00740 
00741 /********************************************************/
00742 /*                                                      */
00743 /*                     moveDCToCenter                   */
00744 /*                                                      */
00745 /********************************************************/
00746 
00747 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
00748           in the image center.
00749 
00750     FFTW produces fourier images where the DC component (origin of
00751     fourier space) is located in the upper left corner of the
00752     image. The quadrants are placed like this (using a 4x4 image for
00753     example):
00754 
00755     \code
00756             DC 4 3 3
00757              4 4 3 3
00758              1 1 2 2
00759              1 1 2 2
00760     \endcode
00761 
00762     After applying the function, the quadrants are at their usual places:
00763 
00764     \code
00765             2 2  1 1
00766             2 2  1 1
00767             3 3 DC 4
00768             3 3  4 4
00769     \endcode
00770 
00771     This transformation can be reversed by \ref moveDCToUpperLeft().
00772     Note that the transformation must not be executed in place - input
00773     and output images must be different.
00774 
00775     <b> Declarations:</b>
00776 
00777     pass arguments explicitly:
00778     \code
00779     namespace vigra {
00780         template <class SrcImageIterator, class SrcAccessor,
00781           class DestImageIterator, class DestAccessor>
00782         void moveDCToCenter(SrcImageIterator src_upperleft,
00783                                SrcImageIterator src_lowerright, SrcAccessor sa,
00784                                DestImageIterator dest_upperleft, DestAccessor da)
00785     }
00786     \endcode
00787 
00788 
00789     use argument objects in conjuction with \ref ArgumentObjectFactories:
00790     \code
00791     namespace vigra {
00792         template <class SrcImageIterator, class SrcAccessor,
00793                   class DestImageIterator, class DestAccessor>
00794         inline void moveDCToCenter(
00795             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00796             pair<DestImageIterator, DestAccessor> dest)
00797     }
00798     \endcode
00799 
00800     <b> Usage:</b>
00801 
00802         <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00803         Namespace: vigra
00804 
00805     \code
00806     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
00807     ... // fill image with data
00808 
00809     // create a plan for optimal performance
00810     fftwnd_plan forwardPlan=
00811         fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE );
00812 
00813     // calculate FFT
00814     fftwnd_one(forwardPlan, spatial.begin(), fourier.begin());
00815 
00816     vigra::FFTWComplexImage rearrangedFourier(width, height);
00817     moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
00818 
00819     //delete the plan
00820     fftwnd_destroy_plan(forwardPlan);
00821     \endcode
00822 */
00823 template <class SrcImageIterator, class SrcAccessor,
00824           class DestImageIterator, class DestAccessor>
00825 void moveDCToCenter(SrcImageIterator src_upperleft,
00826                                SrcImageIterator src_lowerright, SrcAccessor sa,
00827                                DestImageIterator dest_upperleft, DestAccessor da)
00828 {
00829     int w= src_lowerright.x - src_upperleft.x;
00830     int h= src_lowerright.y - src_upperleft.y;
00831     int w1 = w/2;
00832     int h1 = h/2;
00833     int w2 = (w+1)/2;
00834     int h2 = (h+1)/2;
00835 
00836     // 2. Quadrant  zum 4.
00837     copyImage(srcIterRange(src_upperleft,
00838                            src_upperleft  + Diff2D(w2, h2), sa),
00839               destIter    (dest_upperleft + Diff2D(w1, h1), da));
00840 
00841     // 4. Quadrant zum 2.
00842     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
00843                            src_lowerright, sa),
00844               destIter    (dest_upperleft, da));
00845 
00846     // 1. Quadrant zum 3.
00847     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
00848                            src_upperleft  + Diff2D(w,  h2), sa),
00849               destIter    (dest_upperleft + Diff2D(0,  h1), da));
00850 
00851     // 3. Quadrant zum 1.
00852     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
00853                            src_upperleft  + Diff2D(w2, h), sa),
00854               destIter    (dest_upperleft + Diff2D(w1, 0), da));
00855 }
00856 
00857 template <class SrcImageIterator, class SrcAccessor,
00858           class DestImageIterator, class DestAccessor>
00859 inline void moveDCToCenter(
00860     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00861     pair<DestImageIterator, DestAccessor> dest)
00862 {
00863     moveDCToCenter(src.first, src.second, src.third,
00864                           dest.first, dest.second);
00865 }
00866 
00867 /********************************************************/
00868 /*                                                      */
00869 /*                   moveDCToUpperLeft                  */
00870 /*                                                      */
00871 /********************************************************/
00872 
00873 /** \brief Rearrange the quadrants of a Fourier image so that the origin is
00874           in the image's upper left.
00875 
00876      This function is the inversion of \ref moveDCToCenter(). See there
00877      for declarations and a usage example.
00878 
00879 */
00880 template <class SrcImageIterator, class SrcAccessor,
00881           class DestImageIterator, class DestAccessor>
00882 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
00883                                SrcImageIterator src_lowerright, SrcAccessor sa,
00884                                DestImageIterator dest_upperleft, DestAccessor da)
00885 {
00886     int w= src_lowerright.x - src_upperleft.x;
00887     int h= src_lowerright.y - src_upperleft.y;
00888     int w2 = w/2;
00889     int h2 = h/2;
00890     int w1 = (w+1)/2;
00891     int h1 = (h+1)/2;
00892 
00893     // 2. Quadrant  zum 4.
00894     copyImage(srcIterRange(src_upperleft,
00895                            src_upperleft  + Diff2D(w2, h2), sa),
00896               destIter    (dest_upperleft + Diff2D(w1, h1), da));
00897 
00898     // 4. Quadrant zum 2.
00899     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
00900                            src_lowerright, sa),
00901               destIter    (dest_upperleft, da));
00902 
00903     // 1. Quadrant zum 3.
00904     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
00905                            src_upperleft  + Diff2D(w,  h2), sa),
00906               destIter    (dest_upperleft + Diff2D(0,  h1), da));
00907 
00908     // 3. Quadrant zum 1.
00909     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
00910                            src_upperleft  + Diff2D(w2, h), sa),
00911               destIter    (dest_upperleft + Diff2D(w1, 0), da));
00912 }
00913 
00914 template <class SrcImageIterator, class SrcAccessor,
00915           class DestImageIterator, class DestAccessor>
00916 inline void moveDCToUpperLeft(
00917     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00918     pair<DestImageIterator, DestAccessor> dest)
00919 {
00920     moveDCToUpperLeft(src.first, src.second, src.third,
00921                                           dest.first, dest.second);
00922 }
00923 
00924 /********************************************************/
00925 /*                                                      */
00926 /*                   applyFourierFilter                 */
00927 /*                                                      */
00928 /********************************************************/
00929 
00930 /** \brief Apply a filter (defined in the frequency domain) to an image.
00931 
00932     After transferring the image into the frequency domain, it is
00933     multiplied pixel-wise with the filter and transformed back. The
00934     result is always put into the given FFTWComplexImage
00935     <TT>destImg</TT> which must have the right size. Finally, the
00936     result will be normalized to compensate for the two FFTs.
00937 
00938     The input and filter images can be scalar images (more precisely,
00939     <TT>SrcAccessor::value_type</TT> must be scalar) or
00940     FFTWComplexImages (in this case, one conversion is saved).
00941 
00942     The DC entry of the filter must be in the upper left, which is the
00943     position where FFTW expects it (see \ref moveDCToUpperLeft()).
00944 
00945     You can optionally pass two fftwnd_plans as last parameters, the
00946     forward plan and the in-place backward plan. Both must have been
00947     created for the right image size (see sample code).
00948 
00949     <b> Declarations:</b>
00950 
00951     pass arguments explicitly:
00952     \code
00953     namespace vigra {
00954         template <class SrcImageIterator, class SrcAccessor,
00955             class FilterImageIterator, class FilterAccessor>
00956         void applyFourierFilter(SrcImageIterator srcUpperLeft,
00957                                 SrcImageIterator srcLowerRight, SrcAccessor sa,
00958                                 FilterImageIterator filterUpperLeft, FilterAccessor fa,
00959                                 FFTWComplexImage & destImg)
00960 
00961         template <class SrcImageIterator, class SrcAccessor,
00962             class FilterImageIterator, class FilterAccessor>
00963         void applyFourierFilter(SrcImageIterator srcUpperLeft,
00964                                 SrcImageIterator srcLowerRight, SrcAccessor sa,
00965                                 FilterImageIterator filterUpperLeft, FilterAccessor fa,
00966                                 FFTWComplexImage & destImg,
00967                                 const fftwnd_plan & forwardPlan, const fftwnd_plan & backwardPlan)
00968     }
00969     \endcode
00970 
00971     use argument objects in conjuction with \ref ArgumentObjectFactories:
00972     \code
00973     namespace vigra {
00974         template <class SrcImageIterator, class SrcAccessor,
00975             class FilterImageIterator, class FilterAccessor>
00976         inline
00977         void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00978                                 pair<FilterImageIterator, FilterAccessor> filter,
00979                                 FFTWComplexImage &destImg)
00980 
00981         template <class SrcImageIterator, class SrcAccessor,
00982             class FilterImageIterator, class FilterAccessor>
00983         inline
00984         void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00985                                 pair<FilterImageIterator, FilterAccessor> filter,
00986                                 FFTWComplexImage &destImg,
00987                                 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
00988     }
00989     \endcode
00990 
00991     <b> Usage:</b>
00992 
00993     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
00994     Namespace: vigra
00995 
00996     \code
00997     // create a Gaussian filter in Fourier space
00998     vigra::FImage gaussFilter(w, h), filter(w, h);
00999     for(int y=0; y<h; ++y)
01000         for(int x=0; x<w; ++x)
01001         {
01002             xx = float(x - w / 2) / w;
01003             yy = float(y - h / 2) / h;
01004 
01005             gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
01006         }
01007 
01008     // applyFourierFilter() expects the filter's DC in the upper left
01009     moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
01010 
01011     vigra::FFTWComplexImage result(w, h);
01012 
01013     vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);
01014     \endcode
01015 
01016     For inspection of the result, \ref FFTWMagnitudeAccessor might be
01017     useful.
01018 
01019     As mentioned, you can pass along two FFTW-plans. This is useful to
01020     speed up the application of many filters of the same size. Look into the FFTW
01021     documentation for details. You can create optimized plans like this:
01022 
01023     \code
01024     // note that the height comes before the width here!
01025     fftwnd_plan forwardPlan= fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_MEASURE );
01026     fftwnd_plan backwardPlan= fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_MEASURE | FFTW_IN_PLACE);
01027 
01028     applyFourierFilter(srcImageRange(image), srcImage(filter), result,
01029                        forwardPlan, backwardPlan);
01030     \endcode
01031 
01032     <em>Note</em>: The creation of the plans in this way can take
01033     several seconds - the FFTW library will measure different possible
01034     implementations to decide which one is the fastest for <em>this
01035     specific environment</em>. The result of those measurements is
01036     called "wisdom" in the FFTW slang and there are functions to save
01037     it to disk.
01038 
01039 */
01040 
01041 // applyFourierFilter versions without fftwnd_plans:
01042 template <class SrcImageIterator, class SrcAccessor,
01043           class FilterImageIterator, class FilterAccessor,
01044           class DestImageIterator, class DestAccessor>
01045 inline
01046 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01047                         pair<FilterImageIterator, FilterAccessor> filter,
01048                         pair<DestImageIterator, DestAccessor> dest)
01049 {
01050     applyFourierFilter(src.first, src.second, src.third,
01051                        filter.first, filter.second,
01052                        dest.first, dest.second);
01053 }
01054 
01055 template <class SrcImageIterator, class SrcAccessor,
01056           class FilterImageIterator, class FilterAccessor,
01057           class DestImageIterator, class DestAccessor>
01058 void applyFourierFilter(SrcImageIterator srcUpperLeft,
01059                         SrcImageIterator srcLowerRight, SrcAccessor sa,
01060                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
01061                         DestImageIterator destUpperLeft, DestAccessor da)
01062 {
01063     // copy real input images into a complex one...
01064     int w= srcLowerRight.x - srcUpperLeft.x;
01065     int h= srcLowerRight.y - srcUpperLeft.y;
01066 
01067     FFTWComplexImage workImage(w, h);
01068     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01069               destImage(workImage, FFTWWriteRealAccessor()));
01070 
01071     // ...and call the impl
01072     FFTWComplexImage const & cworkImage = workImage;
01073     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01074                            filterUpperLeft, fa,
01075                            destUpperLeft, da);
01076 }
01077 
01078 typedef FFTWComplexImage::const_traverser FFTWConstTraverser;
01079 
01080 template <class FilterImageIterator, class FilterAccessor,
01081           class DestImageIterator, class DestAccessor>
01082 inline
01083 void applyFourierFilter(
01084     FFTWComplexImage::const_traverser srcUpperLeft,
01085     FFTWComplexImage::const_traverser srcLowerRight,
01086     FFTWComplexImage::ConstAccessor sa,
01087     FilterImageIterator filterUpperLeft, FilterAccessor fa,
01088     DestImageIterator destUpperLeft, DestAccessor da)
01089 {
01090     int w= srcLowerRight.x - srcUpperLeft.x;
01091 
01092     // test for right memory layout (fftw expects a 2*width*height floats array)
01093     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
01094         applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
01095                                filterUpperLeft, fa,
01096                                destUpperLeft, da);
01097     else
01098     {
01099         FFTWComplexImage workImage(w, h);
01100         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01101                   destImage(workImage));
01102 
01103         FFTWComplexImage const & cworkImage = workImage;
01104         applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01105                                filterUpperLeft, fa,
01106                                destUpperLeft, da);
01107     }
01108 }
01109 
01110 template <class FilterImageIterator, class FilterAccessor,
01111           class DestImageIterator, class DestAccessor>
01112 void applyFourierFilterImpl(
01113     FFTWComplexImage::const_traverser srcUpperLeft,
01114     FFTWComplexImage::const_traverser srcLowerRight,
01115     FFTWComplexImage::ConstAccessor sa,
01116     FilterImageIterator filterUpperLeft, FilterAccessor fa,
01117     DestImageIterator destUpperLeft, DestAccessor da)
01118 {
01119     // create plans and call variant with plan parameters
01120     int w= srcLowerRight.x - srcUpperLeft.x;
01121     int h= srcLowerRight.y - srcUpperLeft.y;
01122 
01123     fftwnd_plan forwardPlan=
01124         fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE );
01125     fftwnd_plan backwardPlan=
01126         fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
01127 
01128     applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
01129                            filterUpperLeft, fa,
01130                            destUpperLeft, da,
01131                            forwardPlan, backwardPlan);
01132 
01133     fftwnd_destroy_plan(forwardPlan);
01134     fftwnd_destroy_plan(backwardPlan);
01135 }
01136 
01137 // applyFourierFilter versions with fftwnd_plans:
01138 template <class SrcImageIterator, class SrcAccessor,
01139           class FilterImageIterator, class FilterAccessor,
01140           class DestImageIterator, class DestAccessor>
01141 inline
01142 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01143                         pair<FilterImageIterator, FilterAccessor> filter,
01144                         pair<DestImageIterator, DestAccessor> dest,
01145                         const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01146 {
01147     applyFourierFilter(src.first, src.second, src.third,
01148                        filter.first, filter.second,
01149                        dest.first, dest.second,
01150                        forwardPlan, backwardPlan);
01151 }
01152 
01153 template <class SrcImageIterator, class SrcAccessor,
01154           class FilterImageIterator, class FilterAccessor,
01155           class DestImageIterator, class DestAccessor>
01156 void applyFourierFilter(SrcImageIterator srcUpperLeft,
01157                         SrcImageIterator srcLowerRight, SrcAccessor sa,
01158                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
01159                         DestImageIterator destUpperLeft, DestAccessor da,
01160                         const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01161 {
01162     int w= srcLowerRight.x - srcUpperLeft.x;
01163     int h= srcLowerRight.y - srcUpperLeft.y;
01164 
01165     FFTWComplexImage workImage(w, h);
01166     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01167               destImage(workImage, FFTWWriteRealAccessor()));
01168 
01169     FFTWComplexImage const & cworkImage = workImage;
01170     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01171                            filterUpperLeft, fa,
01172                            destUpperLeft, da,
01173                            forwardPlan, backwardPlan);
01174 }
01175 
01176 template <class FilterImageIterator, class FilterAccessor,
01177           class DestImageIterator, class DestAccessor>
01178 inline
01179 void applyFourierFilter(
01180     FFTWComplexImage::const_traverser srcUpperLeft,
01181     FFTWComplexImage::const_traverser srcLowerRight,
01182     FFTWComplexImage::ConstAccessor sa,
01183     FilterImageIterator filterUpperLeft, FilterAccessor fa,
01184     DestImageIterator destUpperLeft, DestAccessor da,
01185     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01186 {
01187     applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
01188                            filterUpperLeft, fa,
01189                            destUpperLeft, da,
01190                            forwardPlan, backwardPlan);
01191 }
01192 
01193 template <class FilterImageIterator, class FilterAccessor,
01194           class DestImageIterator, class DestAccessor>
01195 void applyFourierFilterImpl(
01196     FFTWComplexImage::const_traverser srcUpperLeft,
01197     FFTWComplexImage::const_traverser srcLowerRight,
01198     FFTWComplexImage::ConstAccessor sa,
01199     FilterImageIterator filterUpperLeft, FilterAccessor fa,
01200     DestImageIterator destUpperLeft, DestAccessor da,
01201     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01202 {
01203     FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
01204 
01205     // FFT from srcImage to complexResultImg
01206     fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)),
01207                complexResultImg.begin());
01208 
01209     // convolve in freq. domain (in complexResultImg)
01210     combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
01211                      destImage(complexResultImg), std::multiplies<FFTWComplex>());
01212 
01213     // FFT back into spatial domain (inplace in complexResultImg)
01214     fftwnd_one(backwardPlan, complexResultImg.begin(), 0);
01215 
01216     typedef typename
01217         NumericTraits<typename DestAccessor::value_type>::isScalar
01218         isScalarResult;
01219 
01220     // normalization (after FFTs), maybe stripping imaginary part
01221     applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
01222                                         isScalarResult());
01223 }
01224 
01225 template <class DestImageIterator, class DestAccessor>
01226 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
01227                                          DestImageIterator destUpperLeft,
01228                                          DestAccessor da,
01229                                          VigraFalseType)
01230 {
01231     double normFactor= 1.0/(srcImage.width() * srcImage.height());
01232 
01233     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
01234     {
01235         DestImageIterator dIt= destUpperLeft;
01236         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
01237         {
01238             da.setComponent(srcImage(x, y).re*normFactor, dIt, 0);
01239             da.setComponent(srcImage(x, y).im*normFactor, dIt, 1);
01240         }
01241     }
01242 }
01243 
01244 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
01245         FFTWComplexImage::traverser destUpperLeft,
01246         FFTWComplexImage::Accessor da,
01247         VigraFalseType)
01248 {
01249     transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
01250                    linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height())));
01251 }
01252 
01253 template <class DestImageIterator, class DestAccessor>
01254 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
01255                                          DestImageIterator destUpperLeft,
01256                                          DestAccessor da,
01257                                          VigraTrueType)
01258 {
01259     double normFactor= 1.0/(srcImage.width() * srcImage.height());
01260 
01261     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
01262     {
01263         DestImageIterator dIt= destUpperLeft;
01264         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
01265             da.set(srcImage(x, y).re*normFactor, dIt);
01266     }
01267 }
01268 
01269 /**********************************************************/
01270 /*                                                        */
01271 /*                applyFourierFilterFamily                */
01272 /*                                                        */
01273 /**********************************************************/
01274 
01275 /** \brief Apply an array of filters (defined in the frequency domain) to an image.
01276 
01277     This provides the same functionality as \ref applyFourierFilter(),
01278     but applying several filters at once allows to avoid
01279     repeated Fourier transforms of the source image.
01280 
01281     Filters and result images must be stored in \ref vigra::ImageArray data
01282     structures. In contrast to \ref applyFourierFilter(), this function adjusts
01283     the size of the result images and the the length of the array.
01284 
01285     <b> Declarations:</b>
01286 
01287     pass arguments explicitly:
01288     \code
01289     namespace vigra {
01290         template <class SrcImageIterator, class SrcAccessor, class FilterType>
01291         void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
01292                                       SrcImageIterator srcLowerRight, SrcAccessor sa,
01293                                       const ImageArray<FilterType> &filters,
01294                                       ImageArray<FFTWComplexImage> &results)
01295 
01296         template <class SrcImageIterator, class SrcAccessor, class FilterType>
01297         void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
01298                                       SrcImageIterator srcLowerRight, SrcAccessor sa,
01299                                       const ImageArray<FilterType> &filters,
01300                                       ImageArray<FFTWComplexImage> &results,
01301                                       const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01302     }
01303     \endcode
01304 
01305     use argument objects in conjuction with \ref ArgumentObjectFactories:
01306     \code
01307     namespace vigra {
01308         template <class SrcImageIterator, class SrcAccessor, class FilterType>
01309         inline
01310         void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01311                                       const ImageArray<FilterType> &filters,
01312                                       ImageArray<FFTWComplexImage> &results)
01313 
01314         template <class SrcImageIterator, class SrcAccessor, class FilterType>
01315         inline
01316         void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01317                                       const ImageArray<FilterType> &filters,
01318                                       ImageArray<FFTWComplexImage> &results,
01319                                       const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01320     }
01321     \endcode
01322 
01323     <b> Usage:</b>
01324 
01325     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
01326     Namespace: vigra
01327 
01328     \code
01329     // assuming the presence of a real-valued image named "image" to
01330     // be filtered in this example
01331 
01332     vigra::ImageArray<vigra::FImage> filters(16, image.size());
01333 
01334     for (int i=0; i<filters.size(); i++)
01335          // create some meaningful filters here
01336          createMyFilterOfScale(i, destImage(filters[i]));
01337 
01338     vigra::ImageArray<vigra::FFTWComplexImage> results();
01339 
01340     vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
01341     \endcode
01342 
01343     For details about the FFTW-plans, see the \ref
01344     applyFourierFilter() example.
01345 */
01346 
01347 // applyFourierFilterFamily versions without fftwnd_plans:
01348 template <class SrcImageIterator, class SrcAccessor,
01349           class FilterType, class DestImage>
01350 inline
01351 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01352                               const ImageArray<FilterType> &filters,
01353                               ImageArray<DestImage> &results)
01354 {
01355     applyFourierFilterFamily(src.first, src.second, src.third,
01356                              filters, results);
01357 }
01358 
01359 template <class SrcImageIterator, class SrcAccessor,
01360           class FilterType, class DestImage>
01361 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
01362                               SrcImageIterator srcLowerRight, SrcAccessor sa,
01363                               const ImageArray<FilterType> &filters,
01364                               ImageArray<DestImage> &results)
01365 {
01366     int w= srcLowerRight.x - srcUpperLeft.x;
01367     int h= srcLowerRight.y - srcUpperLeft.y;
01368 
01369     FFTWComplexImage workImage(w, h);
01370     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01371               destImage(workImage, FFTWWriteRealAccessor()));
01372 
01373     FFTWComplexImage const & cworkImage = workImage;
01374     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01375                                  filters, results);
01376 }
01377 
01378 template <class FilterType, class DestImage>
01379 inline
01380 void applyFourierFilterFamily(
01381     FFTWComplexImage::const_traverser srcUpperLeft,
01382     FFTWComplexImage::const_traverser srcLowerRight,
01383     FFTWComplexImage::ConstAccessor sa,
01384     const ImageArray<FilterType> &filters,
01385     ImageArray<DestImage> &results)
01386 {
01387     applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
01388                                  filters, results);
01389 }
01390 
01391 template <class FilterType, class DestImage>
01392 void applyFourierFilterFamilyImpl(
01393     FFTWComplexImage::const_traverser srcUpperLeft,
01394     FFTWComplexImage::const_traverser srcLowerRight,
01395     FFTWComplexImage::ConstAccessor sa,
01396     const ImageArray<FilterType> &filters,
01397     ImageArray<DestImage> &results)
01398 {
01399     int w= srcLowerRight.x - srcUpperLeft.x;
01400     int h= srcLowerRight.y - srcUpperLeft.y;
01401 
01402     fftwnd_plan forwardPlan=
01403         fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE );
01404     fftwnd_plan backwardPlan=
01405         fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
01406 
01407     applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
01408                                  filters, results,
01409                                  forwardPlan, backwardPlan);
01410 
01411     fftwnd_destroy_plan(forwardPlan);
01412     fftwnd_destroy_plan(backwardPlan);
01413 }
01414 
01415 // applyFourierFilterFamily versions with fftwnd_plans:
01416 template <class SrcImageIterator, class SrcAccessor,
01417           class FilterType, class DestImage>
01418 inline
01419 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
01420                               const ImageArray<FilterType> &filters,
01421                               ImageArray<DestImage> &results,
01422                               const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01423 {
01424     applyFourierFilterFamily(src.first, src.second, src.third,
01425                                  filters, results,
01426                                  forwardPlan, backwardPlan);
01427 }
01428 
01429 template <class SrcImageIterator, class SrcAccessor,
01430           class FilterType, class DestImage>
01431 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
01432                               SrcImageIterator srcLowerRight, SrcAccessor sa,
01433                               const ImageArray<FilterType> &filters,
01434                               ImageArray<DestImage> &results,
01435                               const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01436 {
01437     int w= srcLowerRight.x - srcUpperLeft.x;
01438     int h= srcLowerRight.y - srcUpperLeft.y;
01439 
01440     FFTWComplexImage workImage(w, h);
01441     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01442               destImage(workImage, FFTWWriteRealAccessor()));
01443 
01444     FFTWComplexImage const & cworkImage = workImage;
01445     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01446                                  filters, results,
01447                                  forwardPlan, backwardPlan);
01448 }
01449 
01450 template <class FilterType, class DestImage>
01451 inline
01452 void applyFourierFilterFamily(
01453     FFTWComplexImage::const_traverser srcUpperLeft,
01454     FFTWComplexImage::const_traverser srcLowerRight,
01455     FFTWComplexImage::ConstAccessor sa,
01456     const ImageArray<FilterType> &filters,
01457     ImageArray<DestImage> &results,
01458     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01459 {
01460     int w= srcLowerRight.x - srcUpperLeft.x;
01461 
01462     // test for right memory layout (fftw expects a 2*width*height floats array)
01463     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
01464         applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
01465                                      filters, results,
01466                                      forwardPlan, backwardPlan);
01467     else
01468     {
01469         FFTWComplexImage workImage(w, h);
01470         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
01471                   destImage(workImage));
01472 
01473         FFTWComplexImage const & cworkImage = workImage;
01474         applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
01475                                      filters, results,
01476                                      forwardPlan, backwardPlan);
01477     }
01478 }
01479 
01480 template <class FilterType, class DestImage>
01481 void applyFourierFilterFamilyImpl(
01482     FFTWComplexImage::const_traverser srcUpperLeft,
01483     FFTWComplexImage::const_traverser srcLowerRight,
01484     FFTWComplexImage::ConstAccessor sa,
01485     const ImageArray<FilterType> &filters,
01486     ImageArray<DestImage> &results,
01487     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
01488 {
01489     // make sure the filter images have the right dimensions
01490     vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
01491                        "applyFourierFilterFamily called with src image size != filters.imageSize()!");
01492 
01493     // make sure the result image array has the right dimensions
01494     results.resize(filters.size());
01495     results.resizeImages(filters.imageSize());
01496 
01497     // FFT from srcImage to freqImage
01498     int w= srcLowerRight.x - srcUpperLeft.x;
01499     int h= srcLowerRight.y - srcUpperLeft.y;
01500 
01501     FFTWComplexImage freqImage(w, h);
01502     FFTWComplexImage result(w, h);
01503 
01504     fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)), freqImage.begin());
01505 
01506     typedef typename
01507         NumericTraits<typename DestImage::Accessor::value_type>::isScalar
01508         isScalarResult;
01509 
01510     // convolve with filters in freq. domain
01511     for (unsigned int i= 0;  i < filters.size(); i++)
01512     {
01513         combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
01514                          destImage(result), std::multiplies<FFTWComplex>());
01515 
01516         // FFT back into spatial domain (inplace in destImage)
01517         fftwnd_one(backwardPlan, result.begin(), 0);
01518 
01519         // normalization (after FFTs), maybe stripping imaginary part
01520         applyFourierFilterImplNormalization(result,
01521                                             results[i].upperLeft(), results[i].accessor(),
01522                                             isScalarResult());
01523     }
01524 }
01525 
01526 //@}
01527 
01528 } // namespace vigra
01529 
01530 #endif // VIGRA_FFTW_HXX

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

html generated using doxygen and Python
VIGRA 1.2.0 (7 Aug 2003)