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