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