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

| vigra/fftw.hxx |
--- | --- | ---

    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.3.2, Jan 27 2005 )                                    */
    00008 /*    You may use, modify, and distribute this software according       */
    00009 /*    to the terms stated in the LICENSE file included in               */
    00010 /*    the VIGRA distribution.                                           */
    00011 /*                                                                      */
    00012 /*    The VIGRA Website is                                              */
    00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00014 /*    Please direct questions, bug reports, and contributions to        */
    00015 /*        koethe@informatik.uni-hamburg.de                              */
    00016 /*                                                                      */
    00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
    00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
    00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
    00020 /*                                                                      */
    00021 /************************************************************************/
    00022
    00023 #ifndef VIGRA_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/transformimage.hxx"
    00031 #include "vigra/combineimages.hxx"
    00032 #include "vigra/numerictraits.hxx"
    00033 #include "vigra/imagecontainer.hxx"
    00034 #include <fftw.h>
    00035
    00036 namespace vigra {
    00037
    00038 /********************************************************/
    00039 /*                                                      */
    00040 /*                    FFTWComplex                       */
    00041 /*                                                      */
    00042 /********************************************************/
    00043
    00044 /* documentation: see fftw3.hxx
    00045 */
    00046 class FFTWComplex
    00047 : public fftw_complex
    00048 {
    00049   public:
    00050         /** The complex' component type, as defined in '<TT>fftw.h</TT>'
    00051         */
    00052     typedef fftw_real value_type;
    00053
    00054         /** reference type (result of operator[])
    00055         */
    00056     typedef fftw_real & reference;
    00057
    00058         /** const reference type (result of operator[] const)
    00059         */
    00060     typedef fftw_real const & const_reference;
    00061
    00062         /** const reference type (result of operator[] const)
    00063         */
    00064     typedef fftw_real * iterator;
    00065
    00066         /** const reference type (result of operator[] const)
    00067         */
    00068     typedef fftw_real const * const_iterator;
    00069
    00070         /** Construct from real and imaginary part.
    00071             Default: 0.
    00072         */
    00073     FFTWComplex(value_type const & ire = 0.0, value_type const & iim = 0.0)
    00074     {
    00075         re = ire;
    00076         im = iim;
    00077     }
    00078
    00079         /** Copy constructor.
    00080         */
    00081     FFTWComplex(FFTWComplex const & o)
    00082     : fftw_complex(o)
    00083     {}
    00084
    00085         /** Construct from plain <TT>fftw_complex</TT>.
    00086         */
    00087     FFTWComplex(fftw_complex const & o)
    00088     : fftw_complex(o)
    00089     {}
    00090
    00091         /** Construct from TinyVector.
    00092         */
    00093     template <class T>
    00094     FFTWComplex(TinyVector<T, 2> const & o)
    00095     {
    00096         re = o[0];
    00097         im = o[1];
    00098     }
    00099
    00100         /** Assignment.
    00101         */
    00102     FFTWComplex& operator=(FFTWComplex const & o)
    00103     {
    00104         re = o.re;
    00105         im = o.im;
    00106         return *this;
    00107     }
    00108
    00109         /** Assignment.
    00110         */
    00111     FFTWComplex& operator=(fftw_complex const & o)
    00112     {
    00113         re = o.re;
    00114         im = o.im;
    00115         return *this;
    00116     }
    00117
    00118         /** Assignment.
    00119         */
    00120     FFTWComplex& operator=(fftw_real const & o)
    00121     {
    00122         re = o;
    00123         im = 0.0;
    00124         return *this;
    00125     }
    00126
    00127         /** Assignment.
    00128         */
    00129     template <class T>
    00130     FFTWComplex& operator=(TinyVector<T, 2> const & o)
    00131     {
    00132         re = o[0];
    00133         im = o[1];
    00134         return *this;
    00135     }
    00136
    00137         /** Unary negation.
    00138         */
    00139     FFTWComplex operator-() const
    00140         { return FFTWComplex(-re, -im); }
    00141
    00142         /** Squared magnitude x*conj(x)
    00143         */
    00144     value_type squaredMagnitude() const
    00145         { return c_re(*this)*c_re(*this)+c_im(*this)*c_im(*this); }
    00146
    00147         /** Magnitude (length of radius vector).
    00148         */
    00149     value_type magnitude() const
    00150         { return VIGRA_CSTD::sqrt(squaredMagnitude()); }
    00151
    00152         /** Phase angle.
    00153         */
    00154     value_type phase() const
    00155         { return VIGRA_CSTD::atan2(c_im(*this),c_re(*this)); }
    00156
    00157         /** Access components as if number were a vector.
    00158         */
    00159     reference operator[](int i)
    00160         { return (&re)[i]; }
    00161
    00162         /** Read components as if number were a vector.
    00163         */
    00164     const_reference operator[](int i) const
    00165         { return (&re)[i]; }
    00166
    00167         /** Length of complex number (always 2).
    00168         */
    00169     int size() const
    00170         { return 2; }
    00171
    00172     iterator begin()
    00173         { return &re; }
    00174
    00175     iterator end()
    00176         { return &re + 2; }
    00177
    00178     const_iterator begin() const
    00179         { return &re; }
    00180
    00181     const_iterator end() const
    00182         { return &re + 2; }
    00183 };
    00184
    00185 /********************************************************/
    00186 /*                                                      */
    00187 /*                    FFTWComplex Traits                */
    00188 /*                                                      */
    00189 /********************************************************/
    00190
    00191 /* documentation: see fftw3.hxx
    00192 */
    00193 template<>
    00194 struct NumericTraits<fftw_complex>
    00195 {
    00196     typedef fftw_complex Type;
    00197     typedef fftw_complex Promote;
    00198     typedef fftw_complex RealPromote;
    00199     typedef fftw_complex ComplexPromote;
    00200     typedef fftw_real    ValueType;
    00201     typedef VigraFalseType isIntegral;
    00202     typedef VigraFalseType isScalar;
    00203     typedef VigraFalseType isOrdered;
    00204     typedef VigraTrueType  isComplex;
    00205
    00206     static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); }
    00207     static FFTWComplex one() { return FFTWComplex(1.0, 0.0); }
    00208     static FFTWComplex nonZero() { return one(); }
    00209
    00210     static const Promote & toPromote(const Type & v) { return v; }
    00211     static const RealPromote & toRealPromote(const Type & v) { return v; }
    00212     static const Type & fromPromote(const Promote & v) { return v; }
    00213     static const Type & fromRealPromote(const RealPromote & v) { return v; }
    00214 };
    00215
    00216 template<>
    00217 struct NumericTraits<FFTWComplex>
    00218 : public NumericTraits<fftw_complex>
    00219 {
    00220     typedef FFTWComplex Type;
    00221     typedef FFTWComplex Promote;
    00222     typedef FFTWComplex RealPromote;
    00223     typedef FFTWComplex ComplexPromote;
    00224     typedef fftw_real    ValueType;
    00225     typedef VigraFalseType isIntegral;
    00226     typedef VigraFalseType isScalar;
    00227     typedef VigraFalseType isOrdered;
    00228     typedef VigraTrueType  isComplex;
    00229 };
    00230
    00231 template <>
    00232 struct PromoteTraits<fftw_complex, fftw_complex>
    00233 {
    00234     typedef fftw_complex Promote;
    00235 };
    00236
    00237 template <>
    00238 struct PromoteTraits<fftw_complex, double>
    00239 {
    00240     typedef fftw_complex Promote;
    00241 };
    00242
    00243 template <>
    00244 struct PromoteTraits<double, fftw_complex>
    00245 {
    00246     typedef fftw_complex Promote;
    00247 };
    00248
    00249 template <>
    00250 struct PromoteTraits<FFTWComplex, FFTWComplex>
    00251 {
    00252     typedef FFTWComplex Promote;
    00253 };
    00254
    00255 template <>
    00256 struct PromoteTraits<FFTWComplex, double>
    00257 {
    00258     typedef FFTWComplex Promote;
    00259 };
    00260
    00261 template <>
    00262 struct PromoteTraits<double, FFTWComplex>
    00263 {
    00264     typedef FFTWComplex Promote;
    00265 };
    00266
    00267
    00268 /********************************************************/
    00269 /*                                                      */
    00270 /*                    FFTWComplex Operations            */
    00271 /*                                                      */
    00272 /********************************************************/
    00273
    00274 /* documentation: see fftw3.hxx
    00275 */
    00276 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) {
    00277     return c_re(a) == c_re(b) && c_im(a) == c_im(b);
    00278 }
    00279
    00280 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) {
    00281     return c_re(a) != c_re(b) || c_im(a) != c_im(b);
    00282 }
    00283
    00284 inline FFTWComplex &operator +=(FFTWComplex &a, const FFTWComplex &b) {
    00285     c_re(a) += c_re(b);
    00286     c_im(a) += c_im(b);
    00287     return a;
    00288 }
    00289
    00290 inline FFTWComplex &operator -=(FFTWComplex &a, const FFTWComplex &b) {
    00291     c_re(a) -= c_re(b);
    00292     c_im(a) -= c_im(b);
    00293     return a;
    00294 }
    00295
    00296 inline FFTWComplex &operator *=(FFTWComplex &a, const FFTWComplex &b) {
    00297     FFTWComplex::value_type t = c_re(a)*c_re(b)-c_im(a)*c_im(b);
    00298     c_im(a) = c_re(a)*c_im(b)+c_im(a)*c_re(b);
    00299     c_re(a) = t;
    00300     return a;
    00301 }
    00302
    00303 inline FFTWComplex &operator /=(FFTWComplex &a, const FFTWComplex &b) {
    00304     FFTWComplex::value_type sm = b.squaredMagnitude();
    00305     FFTWComplex::value_type t = (c_re(a)*c_re(b)+c_im(a)*c_im(b))/sm;
    00306     c_im(a) = (c_re(b)*c_im(a)-c_re(a)*c_im(b))/sm;
    00307     c_re(a) = t;
    00308     return a;
    00309 }
    00310
    00311 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) {
    00312     c_re(a) *= b;
    00313     c_im(a) *= b;
    00314     return a;
    00315 }
    00316
    00317 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) {
    00318     c_re(a) /= b;
    00319     c_im(a) /= b;
    00320     return a;
    00321 }
    00322
    00323 inline FFTWComplex operator +(FFTWComplex a, const FFTWComplex &b) {
    00324     a += b;
    00325     return a;
    00326 }
    00327
    00328 inline FFTWComplex operator -(FFTWComplex a, const FFTWComplex &b) {
    00329     a -= b;
    00330     return a;
    00331 }
    00332
    00333 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) {
    00334     a *= b;
    00335     return a;
    00336 }
    00337
    00338 inline FFTWComplex operator *(FFTWComplex a, const double &b) {
    00339     a *= b;
    00340     return a;
    00341 }
    00342
    00343 inline FFTWComplex operator *(const double &a, FFTWComplex b) {
    00344     b *= a;
    00345     return b;
    00346 }
    00347
    00348 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) {
    00349     a /= b;
    00350     return a;
    00351 }
    00352
    00353 inline FFTWComplex operator /(FFTWComplex a, const double &b) {
    00354     a /= b;
    00355     return a;
    00356 }
    00357
    00358 using VIGRA_CSTD::abs;
    00359
    00360 inline FFTWComplex::value_type abs(const FFTWComplex &a)
    00361 {
    00362     return a.magnitude();
    00363 }
    00364
    00365 inline FFTWComplex conj(const FFTWComplex &a)
    00366 {
    00367     return FFTWComplex(a.re, -a.im);
    00368 }
    00369
    00370
    00371
    00372 /********************************************************/
    00373 /*                                                      */
    00374 /*                      FFTWRealImage                   */
    00375 /*                                                      */
    00376 /********************************************************/
    00377
    00378 /* documentation: see fftw3.hxx
    00379 */
    00380 typedef BasicImage<fftw_real> FFTWRealImage;
    00381
    00382 /********************************************************/
    00383 /*                                                      */
    00384 /*                     FFTWComplexImage                 */
    00385 /*                                                      */
    00386 /********************************************************/
    00387
    00388 template<>
    00389 struct IteratorTraits<
    00390         BasicImageIterator<FFTWComplex, FFTWComplex **> >
    00391 {
    00392     typedef BasicImageIterator<FFTWComplex, FFTWComplex **>  Iterator;
    00393     typedef Iterator                             iterator;
    00394     typedef iterator::iterator_category          iterator_category;
    00395     typedef iterator::value_type                 value_type;
    00396     typedef iterator::reference                  reference;
    00397     typedef iterator::index_reference            index_reference;
    00398     typedef iterator::pointer                    pointer;
    00399     typedef iterator::difference_type            difference_type;
    00400     typedef iterator::row_iterator               row_iterator;
    00401     typedef iterator::column_iterator            column_iterator;
    00402     typedef VectorAccessor<FFTWComplex>          default_accessor;
    00403     typedef VectorAccessor<FFTWComplex>          DefaultAccessor;
    00404 };
    00405
    00406 template<>
    00407 struct IteratorTraits<
    00408         ConstBasicImageIterator<FFTWComplex, FFTWComplex **> >
    00409 {
    00410     typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **>    Iterator;
    00411     typedef Iterator                             iterator;
    00412     typedef iterator::iterator_category          iterator_category;
    00413     typedef iterator::value_type                 value_type;
    00414     typedef iterator::reference                  reference;
    00415     typedef iterator::index_reference            index_reference;
    00416     typedef iterator::pointer                    pointer;
    00417     typedef iterator::difference_type            difference_type;
    00418     typedef iterator::row_iterator               row_iterator;
    00419     typedef iterator::column_iterator            column_iterator;
    00420     typedef VectorAccessor<FFTWComplex>          default_accessor;
    00421     typedef VectorAccessor<FFTWComplex>          DefaultAccessor;
    00422 };
    00423
    00424 /* documentation: see fftw3.hxx
    00425 */
    00426 typedef BasicImage<FFTWComplex> FFTWComplexImage;
    00427
    00428 /********************************************************/
    00429 /*                                                      */
    00430 /*                  FFTWComplex-Accessors               */
    00431 /*                                                      */
    00432 /********************************************************/
    00433
    00434 /* documentation: see fftw3.hxx
    00435 */
    00436 class FFTWRealAccessor
    00437 {
    00438   public:
    00439
    00440         /// The accessor's value type.
    00441     typedef fftw_real value_type;
    00442
    00443         /// Read real part at iterator position.
    00444     template <class ITERATOR>
    00445     value_type operator()(ITERATOR const & i) const {
    00446         return c_re(*i);
    00447     }
    00448
    00449         /// Read real part at offset from iterator position.
    00450     template <class ITERATOR, class DIFFERENCE>
    00451     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00452         return c_re(i[d]);
    00453     }
    00454
    00455         /// Write real part at iterator position.
    00456     template <class ITERATOR>
    00457     void set(value_type const & v, ITERATOR const & i) const {
    00458         c_re(*i)= v;
    00459     }
    00460
    00461         /// Write real part at offset from iterator position.
    00462     template <class ITERATOR, class DIFFERENCE>
    00463     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
    00464         c_re(i[d])= v;
    00465     }
    00466 };
    00467
    00468 /* documentation: see fftw3.hxx
    00469 */
    00470 class FFTWImaginaryAccessor
    00471 {
    00472   public:
    00473         /// The accessor's value type.
    00474     typedef fftw_real value_type;
    00475
    00476         /// Read imaginary part at iterator position.
    00477     template <class ITERATOR>
    00478     value_type operator()(ITERATOR const & i) const {
    00479         return c_im(*i);
    00480     }
    00481
    00482         /// Read imaginary part at offset from iterator position.
    00483     template <class ITERATOR, class DIFFERENCE>
    00484     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00485         return c_im(i[d]);
    00486     }
    00487
    00488         /// Write imaginary part at iterator position.
    00489     template <class ITERATOR>
    00490     void set(value_type const & v, ITERATOR const & i) const {
    00491         c_im(*i)= v;
    00492     }
    00493
    00494         /// Write imaginary part at offset from iterator position.
    00495     template <class ITERATOR, class DIFFERENCE>
    00496     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
    00497         c_im(i[d])= v;
    00498     }
    00499 };
    00500
    00501 /* documentation: see fftw3.hxx
    00502 */
    00503 class FFTWWriteRealAccessor: public FFTWRealAccessor
    00504 {
    00505   public:
    00506         /// The accessor's value type.
    00507     typedef fftw_real value_type;
    00508
    00509         /** Write real number at iterator position. Set imaginary part
    00510             to 0.
    00511         */
    00512     template <class ITERATOR>
    00513     void set(value_type const & v, ITERATOR const & i) const {
    00514         c_re(*i)= v;
    00515         c_im(*i)= 0;
    00516     }
    00517
    00518         /** Write real number at offset from iterator position. Set imaginary part
    00519             to 0.
    00520         */
    00521     template <class ITERATOR, class DIFFERENCE>
    00522     void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const {
    00523         c_re(i[d])= v;
    00524         c_im(i[d])= 0;
    00525     }
    00526 };
    00527
    00528 /* documentation: see fftw3.hxx
    00529 */
    00530 class FFTWMagnitudeAccessor
    00531 {
    00532   public:
    00533         /// The accessor's value type.
    00534     typedef fftw_real value_type;
    00535
    00536         /// Read magnitude at iterator position.
    00537     template <class ITERATOR>
    00538     value_type operator()(ITERATOR const & i) const {
    00539         return (*i).magnitude();
    00540     }
    00541
    00542         /// Read magnitude at offset from iterator position.
    00543     template <class ITERATOR, class DIFFERENCE>
    00544     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00545         return (i[d]).magnitude();
    00546     }
    00547 };
    00548
    00549 /* documentation: see fftw3.hxx
    00550 */
    00551 class FFTWPhaseAccessor
    00552 {
    00553   public:
    00554         /// The accessor's value type.
    00555     typedef fftw_real value_type;
    00556
    00557         /// Read phase at iterator position.
    00558     template <class ITERATOR>
    00559     value_type operator()(ITERATOR const & i) const {
    00560         return (*i).phase();
    00561     }
    00562
    00563         /// Read phase at offset from iterator position.
    00564     template <class ITERATOR, class DIFFERENCE>
    00565     value_type operator()(ITERATOR const & i, DIFFERENCE d) const {
    00566         return (i[d]).phase();
    00567     }
    00568 };
    00569
    00570 /********************************************************/
    00571 /*                                                      */
    00572 /*                    Fourier Transform                 */
    00573 /*                                                      */
    00574 /********************************************************/
    00575
    00576 /** \page FourierTransformFFTW2 Fast Fourier Transform
    00577
    00578     This documentation describes the deprecated VIGRA interface to
    00579     FFTW 2. Use the \link FourierTransform interface to the newer
    00580     version FFTW 3\endlink instead.
    00581
    00582     VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier
    00583     Transform</a> package to perform Fourier transformations. VIGRA
    00584     provides a wrapper for FFTW's complex number type (FFTWComplex),
    00585     but FFTW's functions are used verbatim. If the image is stored as
    00586     a FFTWComplexImage, a FFT is performed like this:
    00587
    00588     \code
    00589     vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
    00590     ... // fill image with data
    00591
    00592     // create a plan for optimal performance
    00593     fftwnd_plan forwardPlan=
    00594         fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE );
    00595
    00596     // calculate FFT
    00597     fftwnd_one(forwardPlan, spatial.begin(), fourier.begin());
    00598     \endcode
    00599
    00600     Note that in the creation of a plan, the height must be given
    00601     first. Note also that <TT>spatial.begin()</TT> may only be passed
    00602     to <TT>fftwnd_one</TT> if the transform shall be applied to the
    00603     entire image. When you want to retrict operation to an ROI, you
    00604     create a copy of the ROI in an image of appropriate size.
    00605
    00606     More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>.
    00607
    00608     FFTW produces fourier images that have the DC component (the
    00609     origin of the Fourier space) in the upper left corner. Often, one
    00610     wants the origin in the center of the image, so that frequencies
    00611     always increase towards the border of the image. This can be
    00612     achieved by calling \ref moveDCToCenter(). The inverse
    00613     transformation is done by \ref moveDCToUpperLeft().
    00614
    00615     <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br>
    00616     Namespace: vigra
    00617 */
    00618
    00619 /********************************************************/
    00620 /*                                                      */
    00621 /*                     moveDCToCenter                   */
    00622 /*                                                      */
    00623 /********************************************************/
    00624
    00625 /* documentation: see fftw3.hxx
    00626 */
    00627 template <class SrcImageIterator, class SrcAccessor,
    00628           class DestImageIterator, class DestAccessor>
    00629 void moveDCToCenter(SrcImageIterator src_upperleft,
    00630                                SrcImageIterator src_lowerright, SrcAccessor sa,
    00631                                DestImageIterator dest_upperleft, DestAccessor da)
    00632 {
    00633     int w= src_lowerright.x - src_upperleft.x;
    00634     int h= src_lowerright.y - src_upperleft.y;
    00635     int w1 = w/2;
    00636     int h1 = h/2;
    00637     int w2 = (w+1)/2;
    00638     int h2 = (h+1)/2;
    00639
    00640     // 2. Quadrant  zum 4.
    00641     copyImage(srcIterRange(src_upperleft,
    00642                            src_upperleft  + Diff2D(w2, h2), sa),
    00643               destIter    (dest_upperleft + Diff2D(w1, h1), da));
    00644
    00645     // 4. Quadrant zum 2.
    00646     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
    00647                            src_lowerright, sa),
    00648               destIter    (dest_upperleft, da));
    00649
    00650     // 1. Quadrant zum 3.
    00651     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
    00652                            src_upperleft  + Diff2D(w,  h2), sa),
    00653               destIter    (dest_upperleft + Diff2D(0,  h1), da));
    00654
    00655     // 3. Quadrant zum 1.
    00656     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
    00657                            src_upperleft  + Diff2D(w2, h), sa),
    00658               destIter    (dest_upperleft + Diff2D(w1, 0), da));
    00659 }
    00660
    00661 template <class SrcImageIterator, class SrcAccessor,
    00662           class DestImageIterator, class DestAccessor>
    00663 inline void moveDCToCenter(
    00664     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00665     pair<DestImageIterator, DestAccessor> dest)
    00666 {
    00667     moveDCToCenter(src.first, src.second, src.third,
    00668                           dest.first, dest.second);
    00669 }
    00670
    00671 /********************************************************/
    00672 /*                                                      */
    00673 /*                   moveDCToUpperLeft                  */
    00674 /*                                                      */
    00675 /********************************************************/
    00676
    00677 /* documentation: see fftw3.hxx
    00678 */
    00679 template <class SrcImageIterator, class SrcAccessor,
    00680           class DestImageIterator, class DestAccessor>
    00681 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
    00682                                SrcImageIterator src_lowerright, SrcAccessor sa,
    00683                                DestImageIterator dest_upperleft, DestAccessor da)
    00684 {
    00685     int w= src_lowerright.x - src_upperleft.x;
    00686     int h= src_lowerright.y - src_upperleft.y;
    00687     int w2 = w/2;
    00688     int h2 = h/2;
    00689     int w1 = (w+1)/2;
    00690     int h1 = (h+1)/2;
    00691
    00692     // 2. Quadrant  zum 4.
    00693     copyImage(srcIterRange(src_upperleft,
    00694                            src_upperleft  + Diff2D(w2, h2), sa),
    00695               destIter    (dest_upperleft + Diff2D(w1, h1), da));
    00696
    00697     // 4. Quadrant zum 2.
    00698     copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
    00699                            src_lowerright, sa),
    00700               destIter    (dest_upperleft, da));
    00701
    00702     // 1. Quadrant zum 3.
    00703     copyImage(srcIterRange(src_upperleft  + Diff2D(w2, 0),
    00704                            src_upperleft  + Diff2D(w,  h2), sa),
    00705               destIter    (dest_upperleft + Diff2D(0,  h1), da));
    00706
    00707     // 3. Quadrant zum 1.
    00708     copyImage(srcIterRange(src_upperleft  + Diff2D(0,  h2),
    00709                            src_upperleft  + Diff2D(w2, h), sa),
    00710               destIter    (dest_upperleft + Diff2D(w1, 0), da));
    00711 }
    00712
    00713 template <class SrcImageIterator, class SrcAccessor,
    00714           class DestImageIterator, class DestAccessor>
    00715 inline void moveDCToUpperLeft(
    00716     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00717     pair<DestImageIterator, DestAccessor> dest)
    00718 {
    00719     moveDCToUpperLeft(src.first, src.second, src.third,
    00720                                           dest.first, dest.second);
    00721 }
    00722
    00723 /********************************************************/
    00724 /*                                                      */
    00725 /*                   applyFourierFilter                 */
    00726 /*                                                      */
    00727 /********************************************************/
    00728
    00729 /* documentation: see fftw3.hxx
    00730 */
    00731
    00732 // applyFourierFilter versions without fftwnd_plans:
    00733 template <class SrcImageIterator, class SrcAccessor,
    00734           class FilterImageIterator, class FilterAccessor,
    00735           class DestImageIterator, class DestAccessor>
    00736 inline
    00737 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00738                         pair<FilterImageIterator, FilterAccessor> filter,
    00739                         pair<DestImageIterator, DestAccessor> dest)
    00740 {
    00741     applyFourierFilter(src.first, src.second, src.third,
    00742                        filter.first, filter.second,
    00743                        dest.first, dest.second);
    00744 }
    00745
    00746 template <class SrcImageIterator, class SrcAccessor,
    00747           class FilterImageIterator, class FilterAccessor,
    00748           class DestImageIterator, class DestAccessor>
    00749 void applyFourierFilter(SrcImageIterator srcUpperLeft,
    00750                         SrcImageIterator srcLowerRight, SrcAccessor sa,
    00751                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
    00752                         DestImageIterator destUpperLeft, DestAccessor da)
    00753 {
    00754     // copy real input images into a complex one...
    00755     int w= srcLowerRight.x - srcUpperLeft.x;
    00756     int h= srcLowerRight.y - srcUpperLeft.y;
    00757
    00758     FFTWComplexImage workImage(w, h);
    00759     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    00760               destImage(workImage, FFTWWriteRealAccessor()));
    00761
    00762     // ...and call the impl
    00763     FFTWComplexImage const & cworkImage = workImage;
    00764     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    00765                            filterUpperLeft, fa,
    00766                            destUpperLeft, da);
    00767 }
    00768
    00769 typedef FFTWComplexImage::const_traverser FFTWConstTraverser;
    00770
    00771 template <class FilterImageIterator, class FilterAccessor,
    00772           class DestImageIterator, class DestAccessor>
    00773 inline
    00774 void applyFourierFilter(
    00775     FFTWComplexImage::const_traverser srcUpperLeft,
    00776     FFTWComplexImage::const_traverser srcLowerRight,
    00777     FFTWComplexImage::ConstAccessor sa,
    00778     FilterImageIterator filterUpperLeft, FilterAccessor fa,
    00779     DestImageIterator destUpperLeft, DestAccessor da)
    00780 {
    00781     int w= srcLowerRight.x - srcUpperLeft.x;
    00782     int h= srcLowerRight.y - srcUpperLeft.y;
    00783
    00784     // test for right memory layout (fftw expects a 2*width*height floats array)
    00785     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
    00786         applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
    00787                                filterUpperLeft, fa,
    00788                                destUpperLeft, da);
    00789     else
    00790     {
    00791         FFTWComplexImage workImage(w, h);
    00792         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    00793                   destImage(workImage));
    00794
    00795         FFTWComplexImage const & cworkImage = workImage;
    00796         applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    00797                                filterUpperLeft, fa,
    00798                                destUpperLeft, da);
    00799     }
    00800 }
    00801
    00802 template <class FilterImageIterator, class FilterAccessor,
    00803           class DestImageIterator, class DestAccessor>
    00804 void applyFourierFilterImpl(
    00805     FFTWComplexImage::const_traverser srcUpperLeft,
    00806     FFTWComplexImage::const_traverser srcLowerRight,
    00807     FFTWComplexImage::ConstAccessor sa,
    00808     FilterImageIterator filterUpperLeft, FilterAccessor fa,
    00809     DestImageIterator destUpperLeft, DestAccessor da)
    00810 {
    00811     // create plans and call variant with plan parameters
    00812     int w= srcLowerRight.x - srcUpperLeft.x;
    00813     int h= srcLowerRight.y - srcUpperLeft.y;
    00814
    00815     fftwnd_plan forwardPlan=
    00816         fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE );
    00817     fftwnd_plan backwardPlan=
    00818         fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
    00819
    00820     applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
    00821                            filterUpperLeft, fa,
    00822                            destUpperLeft, da,
    00823                            forwardPlan, backwardPlan);
    00824
    00825     fftwnd_destroy_plan(forwardPlan);
    00826     fftwnd_destroy_plan(backwardPlan);
    00827 }
    00828
    00829 // applyFourierFilter versions with fftwnd_plans:
    00830 template <class SrcImageIterator, class SrcAccessor,
    00831           class FilterImageIterator, class FilterAccessor,
    00832           class DestImageIterator, class DestAccessor>
    00833 inline
    00834 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00835                         pair<FilterImageIterator, FilterAccessor> filter,
    00836                         pair<DestImageIterator, DestAccessor> dest,
    00837                         const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    00838 {
    00839     applyFourierFilter(src.first, src.second, src.third,
    00840                        filter.first, filter.second,
    00841                        dest.first, dest.second,
    00842                        forwardPlan, backwardPlan);
    00843 }
    00844
    00845 template <class SrcImageIterator, class SrcAccessor,
    00846           class FilterImageIterator, class FilterAccessor,
    00847           class DestImageIterator, class DestAccessor>
    00848 void applyFourierFilter(SrcImageIterator srcUpperLeft,
    00849                         SrcImageIterator srcLowerRight, SrcAccessor sa,
    00850                         FilterImageIterator filterUpperLeft, FilterAccessor fa,
    00851                         DestImageIterator destUpperLeft, DestAccessor da,
    00852                         const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    00853 {
    00854     int w= srcLowerRight.x - srcUpperLeft.x;
    00855     int h= srcLowerRight.y - srcUpperLeft.y;
    00856
    00857     FFTWComplexImage workImage(w, h);
    00858     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    00859               destImage(workImage, FFTWWriteRealAccessor()));
    00860
    00861     FFTWComplexImage const & cworkImage = workImage;
    00862     applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    00863                            filterUpperLeft, fa,
    00864                            destUpperLeft, da,
    00865                            forwardPlan, backwardPlan);
    00866 }
    00867
    00868 template <class FilterImageIterator, class FilterAccessor,
    00869           class DestImageIterator, class DestAccessor>
    00870 inline
    00871 void applyFourierFilter(
    00872     FFTWComplexImage::const_traverser srcUpperLeft,
    00873     FFTWComplexImage::const_traverser srcLowerRight,
    00874     FFTWComplexImage::ConstAccessor sa,
    00875     FilterImageIterator filterUpperLeft, FilterAccessor fa,
    00876     DestImageIterator destUpperLeft, DestAccessor da,
    00877     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    00878 {
    00879     applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
    00880                            filterUpperLeft, fa,
    00881                            destUpperLeft, da,
    00882                            forwardPlan, backwardPlan);
    00883 }
    00884
    00885 template <class FilterImageIterator, class FilterAccessor,
    00886           class DestImageIterator, class DestAccessor>
    00887 void applyFourierFilterImpl(
    00888     FFTWComplexImage::const_traverser srcUpperLeft,
    00889     FFTWComplexImage::const_traverser srcLowerRight,
    00890     FFTWComplexImage::ConstAccessor sa,
    00891     FilterImageIterator filterUpperLeft, FilterAccessor fa,
    00892     DestImageIterator destUpperLeft, DestAccessor da,
    00893     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    00894 {
    00895     FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
    00896
    00897     // FFT from srcImage to complexResultImg
    00898     fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)),
    00899                complexResultImg.begin());
    00900
    00901     // convolve in freq. domain (in complexResultImg)
    00902     combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
    00903                      destImage(complexResultImg), std::multiplies<FFTWComplex>());
    00904
    00905     // FFT back into spatial domain (inplace in complexResultImg)
    00906     fftwnd_one(backwardPlan, complexResultImg.begin(), 0);
    00907
    00908     typedef typename
    00909         NumericTraits<typename DestAccessor::value_type>::isScalar
    00910         isScalarResult;
    00911
    00912     // normalization (after FFTs), maybe stripping imaginary part
    00913     applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
    00914                                         isScalarResult());
    00915 }
    00916
    00917 template <class DestImageIterator, class DestAccessor>
    00918 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage,
    00919                                          DestImageIterator destUpperLeft,
    00920                                          DestAccessor da,
    00921                                          VigraFalseType)
    00922 {
    00923     double normFactor= 1.0/(srcImage.width() * srcImage.height());
    00924
    00925     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
    00926     {
    00927         DestImageIterator dIt= destUpperLeft;
    00928         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
    00929         {
    00930             da.setComponent(srcImage(x, y).re*normFactor, dIt, 0);
    00931             da.setComponent(srcImage(x, y).im*normFactor, dIt, 1);
    00932         }
    00933     }
    00934 }
    00935
    00936 inline
    00937 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
    00938         FFTWComplexImage::traverser destUpperLeft,
    00939         FFTWComplexImage::Accessor da,
    00940         VigraFalseType)
    00941 {
    00942     transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
    00943                    linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height())));
    00944 }
    00945
    00946 template <class DestImageIterator, class DestAccessor>
    00947 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage,
    00948                                          DestImageIterator destUpperLeft,
    00949                                          DestAccessor da,
    00950                                          VigraTrueType)
    00951 {
    00952     double normFactor= 1.0/(srcImage.width() * srcImage.height());
    00953
    00954     for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
    00955     {
    00956         DestImageIterator dIt= destUpperLeft;
    00957         for(int x= 0; x< srcImage.width(); x++, dIt.x++)
    00958             da.set(srcImage(x, y).re*normFactor, dIt);
    00959     }
    00960 }
    00961
    00962 /**********************************************************/
    00963 /*                                                        */
    00964 /*                applyFourierFilterFamily                */
    00965 /*                                                        */
    00966 /**********************************************************/
    00967
    00968 /* documentation: see fftw3.hxx
    00969 */
    00970
    00971 // applyFourierFilterFamily versions without fftwnd_plans:
    00972 template <class SrcImageIterator, class SrcAccessor,
    00973           class FilterType, class DestImage>
    00974 inline
    00975 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    00976                               const ImageArray<FilterType> &filters,
    00977                               ImageArray<DestImage> &results)
    00978 {
    00979     applyFourierFilterFamily(src.first, src.second, src.third,
    00980                              filters, results);
    00981 }
    00982
    00983 template <class SrcImageIterator, class SrcAccessor,
    00984           class FilterType, class DestImage>
    00985 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
    00986                               SrcImageIterator srcLowerRight, SrcAccessor sa,
    00987                               const ImageArray<FilterType> &filters,
    00988                               ImageArray<DestImage> &results)
    00989 {
    00990     int w= srcLowerRight.x - srcUpperLeft.x;
    00991     int h= srcLowerRight.y - srcUpperLeft.y;
    00992
    00993     FFTWComplexImage workImage(w, h);
    00994     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    00995               destImage(workImage, FFTWWriteRealAccessor()));
    00996
    00997     FFTWComplexImage const & cworkImage = workImage;
    00998     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    00999                                  filters, results);
    01000 }
    01001
    01002 template <class FilterType, class DestImage>
    01003 inline
    01004 void applyFourierFilterFamily(
    01005     FFTWComplexImage::const_traverser srcUpperLeft,
    01006     FFTWComplexImage::const_traverser srcLowerRight,
    01007     FFTWComplexImage::ConstAccessor sa,
    01008     const ImageArray<FilterType> &filters,
    01009     ImageArray<DestImage> &results)
    01010 {
    01011     applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
    01012                                  filters, results);
    01013 }
    01014
    01015 template <class FilterType, class DestImage>
    01016 void applyFourierFilterFamilyImpl(
    01017     FFTWComplexImage::const_traverser srcUpperLeft,
    01018     FFTWComplexImage::const_traverser srcLowerRight,
    01019     FFTWComplexImage::ConstAccessor sa,
    01020     const ImageArray<FilterType> &filters,
    01021     ImageArray<DestImage> &results)
    01022 {
    01023     int w= srcLowerRight.x - srcUpperLeft.x;
    01024     int h= srcLowerRight.y - srcUpperLeft.y;
    01025
    01026     fftwnd_plan forwardPlan=
    01027         fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE );
    01028     fftwnd_plan backwardPlan=
    01029         fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);
    01030
    01031     applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
    01032                                  filters, results,
    01033                                  forwardPlan, backwardPlan);
    01034
    01035     fftwnd_destroy_plan(forwardPlan);
    01036     fftwnd_destroy_plan(backwardPlan);
    01037 }
    01038
    01039 // applyFourierFilterFamily versions with fftwnd_plans:
    01040 template <class SrcImageIterator, class SrcAccessor,
    01041           class FilterType, class DestImage>
    01042 inline
    01043 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
    01044                               const ImageArray<FilterType> &filters,
    01045                               ImageArray<DestImage> &results,
    01046                               const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    01047 {
    01048     applyFourierFilterFamily(src.first, src.second, src.third,
    01049                                  filters, results,
    01050                                  forwardPlan, backwardPlan);
    01051 }
    01052
    01053 template <class SrcImageIterator, class SrcAccessor,
    01054           class FilterType, class DestImage>
    01055 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
    01056                               SrcImageIterator srcLowerRight, SrcAccessor sa,
    01057                               const ImageArray<FilterType> &filters,
    01058                               ImageArray<DestImage> &results,
    01059                               const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    01060 {
    01061     int w= srcLowerRight.x - srcUpperLeft.x;
    01062     int h= srcLowerRight.y - srcUpperLeft.y;
    01063
    01064     FFTWComplexImage workImage(w, h);
    01065     copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    01066               destImage(workImage, FFTWWriteRealAccessor()));
    01067
    01068     FFTWComplexImage const & cworkImage = workImage;
    01069     applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    01070                                  filters, results,
    01071                                  forwardPlan, backwardPlan);
    01072 }
    01073
    01074 template <class FilterType, class DestImage>
    01075 inline
    01076 void applyFourierFilterFamily(
    01077     FFTWComplexImage::const_traverser srcUpperLeft,
    01078     FFTWComplexImage::const_traverser srcLowerRight,
    01079     FFTWComplexImage::ConstAccessor sa,
    01080     const ImageArray<FilterType> &filters,
    01081     ImageArray<DestImage> &results,
    01082     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    01083 {
    01084     int w= srcLowerRight.x - srcUpperLeft.x;
    01085     int h= srcLowerRight.y - srcUpperLeft.y;
    01086
    01087     // test for right memory layout (fftw expects a 2*width*height floats array)
    01088     if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
    01089         applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
    01090                                      filters, results,
    01091                                      forwardPlan, backwardPlan);
    01092     else
    01093     {
    01094         FFTWComplexImage workImage(w, h);
    01095         copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
    01096                   destImage(workImage));
    01097
    01098         FFTWComplexImage const & cworkImage = workImage;
    01099         applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
    01100                                      filters, results,
    01101                                      forwardPlan, backwardPlan);
    01102     }
    01103 }
    01104
    01105 template <class FilterType, class DestImage>
    01106 void applyFourierFilterFamilyImpl(
    01107     FFTWComplexImage::const_traverser srcUpperLeft,
    01108     FFTWComplexImage::const_traverser srcLowerRight,
    01109     FFTWComplexImage::ConstAccessor sa,
    01110     const ImageArray<FilterType> &filters,
    01111     ImageArray<DestImage> &results,
    01112     const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan)
    01113 {
    01114     // make sure the filter images have the right dimensions
    01115     vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
    01116                        "applyFourierFilterFamily called with src image size != filters.imageSize()!");
    01117
    01118     // make sure the result image array has the right dimensions
    01119     results.resize(filters.size());
    01120     results.resizeImages(filters.imageSize());
    01121
    01122     // FFT from srcImage to freqImage
    01123     int w= srcLowerRight.x - srcUpperLeft.x;
    01124     int h= srcLowerRight.y - srcUpperLeft.y;
    01125
    01126     FFTWComplexImage freqImage(w, h);
    01127     FFTWComplexImage result(w, h);
    01128
    01129     fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)), freqImage.begin());
    01130
    01131     typedef typename
    01132         NumericTraits<typename DestImage::Accessor::value_type>::isScalar
    01133         isScalarResult;
    01134
    01135     // convolve with filters in freq. domain
    01136     for (unsigned int i= 0;  i < filters.size(); i++)
    01137     {
    01138         combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]),
    01139                          destImage(result), std::multiplies<FFTWComplex>());
    01140
    01141         // FFT back into spatial domain (inplace in destImage)
    01142         fftwnd_one(backwardPlan, result.begin(), 0);
    01143
    01144         // normalization (after FFTs), maybe stripping imaginary part
    01145         applyFourierFilterImplNormalization(result,
    01146                                             results[i].upperLeft(), results[i].accessor(),
    01147                                             isScalarResult());
    01148     }
    01149 }
    01150
    01151 } // namespace vigra
    01152
    01153 #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.3.2 (27 Jan 2005)