[ 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.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)