[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2004 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_FFTW3_HXX 00037 #define VIGRA_FFTW3_HXX 00038 00039 #include <cmath> 00040 #include <functional> 00041 #include "stdimage.hxx" 00042 #include "copyimage.hxx" 00043 #include "transformimage.hxx" 00044 #include "combineimages.hxx" 00045 #include "numerictraits.hxx" 00046 #include "imagecontainer.hxx" 00047 #include <fftw3.h> 00048 00049 namespace vigra { 00050 00051 typedef double fftw_real; 00052 00053 /********************************************************/ 00054 /* */ 00055 /* FFTWComplex */ 00056 /* */ 00057 /********************************************************/ 00058 00059 /** \brief Wrapper class for the FFTW type '<TT>fftw_complex</TT>'. 00060 00061 This class provides constructors and other member functions 00062 for the C struct '<TT>fftw_complex</TT>'. This struct is the basic 00063 pixel type of the <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a> 00064 library. It inherits the data members '<TT>re</TT>' and '<TT>im</TT>' 00065 that denote the real and imaginary part of the number. In addition it 00066 defines transformations to polar coordinates, 00067 as well as \ref FFTWComplexOperators "arithmetic operators" and 00068 \ref FFTWComplexAccessors "accessors". 00069 00070 FFTWComplex implements the concepts \ref AlgebraicField and 00071 \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt> 00072 and <tt>FFTWComplexImage</tt> are defined. 00073 00074 <b>See also:</b> 00075 <ul> 00076 <li> \ref FFTWComplexTraits 00077 <li> \ref FFTWComplexOperators 00078 <li> \ref FFTWComplexAccessors 00079 </ul> 00080 00081 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00082 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00083 Namespace: vigra 00084 */ 00085 class FFTWComplex 00086 { 00087 fftw_complex data_; 00088 00089 public: 00090 /** The complex' component type, as defined in '<TT>fftw3.h</TT>' 00091 */ 00092 typedef fftw_real value_type; 00093 00094 /** reference type (result of operator[]) 00095 */ 00096 typedef fftw_real & reference; 00097 00098 /** const reference type (result of operator[] const) 00099 */ 00100 typedef fftw_real const & const_reference; 00101 00102 /** iterator type (result of begin() ) 00103 */ 00104 typedef fftw_real * iterator; 00105 00106 /** const iterator type (result of begin() const) 00107 */ 00108 typedef fftw_real const * const_iterator; 00109 00110 /** The norm type (result of magnitde()) 00111 */ 00112 typedef fftw_real NormType; 00113 00114 /** The squared norm type (result of squaredMagnitde()) 00115 */ 00116 typedef fftw_real SquaredNormType; 00117 00118 /** Construct from real and imaginary part. 00119 Default: 0. 00120 */ 00121 FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0) 00122 { 00123 data_[0] = re; 00124 data_[1] = im; 00125 } 00126 00127 /** Copy constructor. 00128 */ 00129 FFTWComplex(FFTWComplex const & o) 00130 { 00131 data_[0] = o.data_[0]; 00132 data_[1] = o.data_[1]; 00133 } 00134 00135 /** Construct from plain <TT>fftw_complex</TT>. 00136 */ 00137 FFTWComplex(fftw_complex const & o) 00138 { 00139 data_[0] = o[0]; 00140 data_[1] = o[1]; 00141 } 00142 00143 /** Construct from TinyVector. 00144 */ 00145 template <class T> 00146 FFTWComplex(TinyVector<T, 2> const & o) 00147 { 00148 data_[0] = o[0]; 00149 data_[1] = o[1]; 00150 } 00151 00152 /** Assignment. 00153 */ 00154 FFTWComplex& operator=(FFTWComplex const & o) 00155 { 00156 data_[0] = o.data_[0]; 00157 data_[1] = o.data_[1]; 00158 return *this; 00159 } 00160 00161 /** Assignment. 00162 */ 00163 FFTWComplex& operator=(fftw_complex const & o) 00164 { 00165 data_[0] = o[0]; 00166 data_[1] = o[1]; 00167 return *this; 00168 } 00169 00170 /** Assignment. 00171 */ 00172 FFTWComplex& operator=(fftw_real const & o) 00173 { 00174 data_[0] = o; 00175 data_[1] = 0.0; 00176 return *this; 00177 } 00178 00179 /** Assignment. 00180 */ 00181 template <class T> 00182 FFTWComplex& operator=(TinyVector<T, 2> const & o) 00183 { 00184 data_[0] = o[0]; 00185 data_[1] = o[1]; 00186 return *this; 00187 } 00188 00189 reference re() 00190 { return data_[0]; } 00191 00192 const_reference re() const 00193 { return data_[0]; } 00194 00195 reference im() 00196 { return data_[1]; } 00197 00198 const_reference im() const 00199 { return data_[1]; } 00200 00201 /** Unary negation. 00202 */ 00203 FFTWComplex operator-() const 00204 { return FFTWComplex(-data_[0], -data_[1]); } 00205 00206 /** Squared magnitude x*conj(x) 00207 */ 00208 SquaredNormType squaredMagnitude() const 00209 { return data_[0]*data_[0]+data_[1]*data_[1]; } 00210 00211 /** Magnitude (length of radius vector). 00212 */ 00213 NormType magnitude() const 00214 { return VIGRA_CSTD::sqrt(squaredMagnitude()); } 00215 00216 /** Phase angle. 00217 */ 00218 value_type phase() const 00219 { return VIGRA_CSTD::atan2(data_[1], data_[0]); } 00220 00221 /** Access components as if number were a vector. 00222 */ 00223 reference operator[](int i) 00224 { return data_[i]; } 00225 00226 /** Read components as if number were a vector. 00227 */ 00228 const_reference operator[](int i) const 00229 { return data_[i]; } 00230 00231 /** Length of complex number (always 2). 00232 */ 00233 int size() const 00234 { return 2; } 00235 00236 iterator begin() 00237 { return data_; } 00238 00239 iterator end() 00240 { return data_ + 2; } 00241 00242 const_iterator begin() const 00243 { return data_; } 00244 00245 const_iterator end() const 00246 { return data_ + 2; } 00247 }; 00248 00249 /********************************************************/ 00250 /* */ 00251 /* FFTWComplexTraits */ 00252 /* */ 00253 /********************************************************/ 00254 00255 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex 00256 00257 The numeric and promote traits for fftw_complex and FFTWComplex follow 00258 the general specifications for \ref NumericPromotionTraits and 00259 \ref AlgebraicField. They are explicitly specialized for the types 00260 involved: 00261 00262 \code 00263 00264 template<> 00265 struct NumericTraits<fftw_complex> 00266 { 00267 typedef fftw_complex Promote; 00268 typedef fftw_complex RealPromote; 00269 typedef fftw_complex ComplexPromote; 00270 typedef fftw_real ValueType; 00271 00272 typedef VigraFalseType isIntegral; 00273 typedef VigraFalseType isScalar; 00274 typedef VigraFalseType isOrdered; 00275 typedef VigraTrueType isComplex; 00276 00277 // etc. 00278 }; 00279 00280 template<> 00281 struct NumericTraits<FFTWComplex> 00282 { 00283 typedef FFTWComplex Promote; 00284 typedef FFTWComplex RealPromote; 00285 typedef FFTWComplex ComplexPromote; 00286 typedef fftw_real ValueType; 00287 00288 typedef VigraFalseType isIntegral; 00289 typedef VigraFalseType isScalar; 00290 typedef VigraFalseType isOrdered; 00291 typedef VigraTrueType isComplex; 00292 00293 // etc. 00294 }; 00295 00296 template<> 00297 struct NormTraits<fftw_complex> 00298 { 00299 typedef fftw_complex Type; 00300 typedef fftw_real SquaredNormType; 00301 typedef fftw_real NormType; 00302 }; 00303 00304 template<> 00305 struct NormTraits<FFTWComplex> 00306 { 00307 typedef FFTWComplex Type; 00308 typedef fftw_real SquaredNormType; 00309 typedef fftw_real NormType; 00310 }; 00311 00312 template <> 00313 struct PromoteTraits<fftw_complex, fftw_complex> 00314 { 00315 typedef fftw_complex Promote; 00316 }; 00317 00318 template <> 00319 struct PromoteTraits<fftw_complex, double> 00320 { 00321 typedef fftw_complex Promote; 00322 }; 00323 00324 template <> 00325 struct PromoteTraits<double, fftw_complex> 00326 { 00327 typedef fftw_complex Promote; 00328 }; 00329 00330 template <> 00331 struct PromoteTraits<FFTWComplex, FFTWComplex> 00332 { 00333 typedef FFTWComplex Promote; 00334 }; 00335 00336 template <> 00337 struct PromoteTraits<FFTWComplex, double> 00338 { 00339 typedef FFTWComplex Promote; 00340 }; 00341 00342 template <> 00343 struct PromoteTraits<double, FFTWComplex> 00344 { 00345 typedef FFTWComplex Promote; 00346 }; 00347 \endcode 00348 00349 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00350 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00351 Namespace: vigra 00352 00353 */ 00354 template<> 00355 struct NumericTraits<fftw_complex> 00356 { 00357 typedef fftw_complex Type; 00358 typedef fftw_complex Promote; 00359 typedef fftw_complex RealPromote; 00360 typedef fftw_complex ComplexPromote; 00361 typedef fftw_real ValueType; 00362 00363 typedef VigraFalseType isIntegral; 00364 typedef VigraFalseType isScalar; 00365 typedef NumericTraits<fftw_real>::isSigned isSigned; 00366 typedef VigraFalseType isOrdered; 00367 typedef VigraTrueType isComplex; 00368 00369 static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); } 00370 static FFTWComplex one() { return FFTWComplex(1.0, 0.0); } 00371 static FFTWComplex nonZero() { return one(); } 00372 00373 static const Promote & toPromote(const Type & v) { return v; } 00374 static const RealPromote & toRealPromote(const Type & v) { return v; } 00375 static const Type & fromPromote(const Promote & v) { return v; } 00376 static const Type & fromRealPromote(const RealPromote & v) { return v; } 00377 }; 00378 00379 template<> 00380 struct NumericTraits<FFTWComplex> 00381 { 00382 typedef FFTWComplex Type; 00383 typedef FFTWComplex Promote; 00384 typedef FFTWComplex RealPromote; 00385 typedef FFTWComplex ComplexPromote; 00386 typedef fftw_real ValueType; 00387 00388 typedef VigraFalseType isIntegral; 00389 typedef VigraFalseType isScalar; 00390 typedef NumericTraits<fftw_real>::isSigned isSigned; 00391 typedef VigraFalseType isOrdered; 00392 typedef VigraTrueType isComplex; 00393 00394 static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); } 00395 static FFTWComplex one() { return FFTWComplex(1.0, 0.0); } 00396 static FFTWComplex nonZero() { return one(); } 00397 00398 static const Promote & toPromote(const Type & v) { return v; } 00399 static const RealPromote & toRealPromote(const Type & v) { return v; } 00400 static const Type & fromPromote(const Promote & v) { return v; } 00401 static const Type & fromRealPromote(const RealPromote & v) { return v; } 00402 }; 00403 00404 template<> 00405 struct NormTraits<fftw_complex> 00406 { 00407 typedef fftw_complex Type; 00408 typedef fftw_real SquaredNormType; 00409 typedef fftw_real NormType; 00410 }; 00411 00412 template<> 00413 struct NormTraits<FFTWComplex> 00414 { 00415 typedef FFTWComplex Type; 00416 typedef fftw_real SquaredNormType; 00417 typedef fftw_real NormType; 00418 }; 00419 00420 template <> 00421 struct PromoteTraits<fftw_complex, fftw_complex> 00422 { 00423 typedef fftw_complex Promote; 00424 }; 00425 00426 template <> 00427 struct PromoteTraits<fftw_complex, double> 00428 { 00429 typedef fftw_complex Promote; 00430 }; 00431 00432 template <> 00433 struct PromoteTraits<double, fftw_complex> 00434 { 00435 typedef fftw_complex Promote; 00436 }; 00437 00438 template <> 00439 struct PromoteTraits<FFTWComplex, FFTWComplex> 00440 { 00441 typedef FFTWComplex Promote; 00442 }; 00443 00444 template <> 00445 struct PromoteTraits<FFTWComplex, double> 00446 { 00447 typedef FFTWComplex Promote; 00448 }; 00449 00450 template <> 00451 struct PromoteTraits<double, FFTWComplex> 00452 { 00453 typedef FFTWComplex Promote; 00454 }; 00455 00456 template<class T> 00457 struct CanSkipInitialization<std::complex<T> > 00458 { 00459 typedef typename CanSkipInitialization<T>::type type; 00460 static const bool value = type::asBool; 00461 }; 00462 00463 template<> 00464 struct CanSkipInitialization<FFTWComplex> 00465 { 00466 typedef CanSkipInitialization<fftw_real>::type type; 00467 static const bool value = type::asBool; 00468 }; 00469 00470 00471 00472 /********************************************************/ 00473 /* */ 00474 /* FFTWComplex Operations */ 00475 /* */ 00476 /********************************************************/ 00477 00478 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex 00479 00480 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00481 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00482 00483 These functions fulfill the requirements of an Algebraic Field. 00484 Return types are determined according to \ref FFTWComplexTraits. 00485 00486 Namespace: vigra 00487 <p> 00488 00489 */ 00490 //@{ 00491 /// equal 00492 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) { 00493 return a.re() == b.re() && a.im() == b.im(); 00494 } 00495 00496 /// not equal 00497 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) { 00498 return a.re() != b.re() || a.im() != b.im(); 00499 } 00500 00501 /// add-assignment 00502 inline FFTWComplex &operator +=(FFTWComplex &a, const FFTWComplex &b) { 00503 a.re() += b.re(); 00504 a.im() += b.im(); 00505 return a; 00506 } 00507 00508 /// subtract-assignment 00509 inline FFTWComplex &operator -=(FFTWComplex &a, const FFTWComplex &b) { 00510 a.re() -= b.re(); 00511 a.im() -= b.im(); 00512 return a; 00513 } 00514 00515 /// multiply-assignment 00516 inline FFTWComplex &operator *=(FFTWComplex &a, const FFTWComplex &b) { 00517 FFTWComplex::value_type t = a.re()*b.re()-a.im()*b.im(); 00518 a.im() = a.re()*b.im()+a.im()*b.re(); 00519 a.re() = t; 00520 return a; 00521 } 00522 00523 /// divide-assignment 00524 inline FFTWComplex &operator /=(FFTWComplex &a, const FFTWComplex &b) { 00525 FFTWComplex::value_type sm = b.squaredMagnitude(); 00526 FFTWComplex::value_type t = (a.re()*b.re()+a.im()*b.im())/sm; 00527 a.im() = (b.re()*a.im()-a.re()*b.im())/sm; 00528 a.re() = t; 00529 return a; 00530 } 00531 00532 /// multiply-assignment with scalar double 00533 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) { 00534 a.re() *= b; 00535 a.im() *= b; 00536 return a; 00537 } 00538 00539 /// divide-assignment with scalar double 00540 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) { 00541 a.re() /= b; 00542 a.im() /= b; 00543 return a; 00544 } 00545 00546 /// addition 00547 inline FFTWComplex operator +(FFTWComplex a, const FFTWComplex &b) { 00548 a += b; 00549 return a; 00550 } 00551 00552 /// subtraction 00553 inline FFTWComplex operator -(FFTWComplex a, const FFTWComplex &b) { 00554 a -= b; 00555 return a; 00556 } 00557 00558 /// multiplication 00559 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) { 00560 a *= b; 00561 return a; 00562 } 00563 00564 /// right multiplication with scalar double 00565 inline FFTWComplex operator *(FFTWComplex a, const double &b) { 00566 a *= b; 00567 return a; 00568 } 00569 00570 /// left multiplication with scalar double 00571 inline FFTWComplex operator *(const double &a, FFTWComplex b) { 00572 b *= a; 00573 return b; 00574 } 00575 00576 /// division 00577 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) { 00578 a /= b; 00579 return a; 00580 } 00581 00582 /// right division with scalar double 00583 inline FFTWComplex operator /(FFTWComplex a, const double &b) { 00584 a /= b; 00585 return a; 00586 } 00587 00588 using VIGRA_CSTD::abs; 00589 00590 /// absolute value (= magnitude) 00591 inline FFTWComplex::value_type abs(const FFTWComplex &a) 00592 { 00593 return a.magnitude(); 00594 } 00595 00596 /// complex conjugate 00597 inline FFTWComplex conj(const FFTWComplex &a) 00598 { 00599 return FFTWComplex(a.re(), -a.im()); 00600 } 00601 00602 /// norm (= magnitude) 00603 inline FFTWComplex::NormType norm(const FFTWComplex &a) 00604 { 00605 return a.magnitude(); 00606 } 00607 00608 /// squared norm (= squared magnitude) 00609 inline FFTWComplex::SquaredNormType squaredNorm(const FFTWComplex &a) 00610 { 00611 return a.squaredMagnitude(); 00612 } 00613 00614 //@} 00615 00616 00617 /** \addtogroup StandardImageTypes 00618 */ 00619 //@{ 00620 00621 /********************************************************/ 00622 /* */ 00623 /* FFTWRealImage */ 00624 /* */ 00625 /********************************************************/ 00626 00627 /** Float (<tt>fftw_real</tt>) image. 00628 00629 The type <tt>fftw_real</tt> is defined as <tt>double</tt> (in FFTW 2 it used to be 00630 either <tt>float</tt> or <tt>double</tt>, as specified during compilation of FFTW). 00631 FFTWRealImage uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and 00632 their const counterparts to access the data. 00633 00634 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00635 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00636 Namespace: vigra 00637 */ 00638 typedef BasicImage<fftw_real> FFTWRealImage; 00639 00640 /********************************************************/ 00641 /* */ 00642 /* FFTWComplexImage */ 00643 /* */ 00644 /********************************************************/ 00645 00646 template<> 00647 struct IteratorTraits< 00648 BasicImageIterator<FFTWComplex, FFTWComplex **> > 00649 { 00650 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> Iterator; 00651 typedef Iterator iterator; 00652 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> mutable_iterator; 00653 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> const_iterator; 00654 typedef iterator::iterator_category iterator_category; 00655 typedef iterator::value_type value_type; 00656 typedef iterator::reference reference; 00657 typedef iterator::index_reference index_reference; 00658 typedef iterator::pointer pointer; 00659 typedef iterator::difference_type difference_type; 00660 typedef iterator::row_iterator row_iterator; 00661 typedef iterator::column_iterator column_iterator; 00662 typedef VectorAccessor<FFTWComplex> default_accessor; 00663 typedef VectorAccessor<FFTWComplex> DefaultAccessor; 00664 typedef VigraTrueType hasConstantStrides; 00665 }; 00666 00667 template<> 00668 struct IteratorTraits< 00669 ConstBasicImageIterator<FFTWComplex, FFTWComplex **> > 00670 { 00671 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> Iterator; 00672 typedef Iterator iterator; 00673 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> mutable_iterator; 00674 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> const_iterator; 00675 typedef iterator::iterator_category iterator_category; 00676 typedef iterator::value_type value_type; 00677 typedef iterator::reference reference; 00678 typedef iterator::index_reference index_reference; 00679 typedef iterator::pointer pointer; 00680 typedef iterator::difference_type difference_type; 00681 typedef iterator::row_iterator row_iterator; 00682 typedef iterator::column_iterator column_iterator; 00683 typedef VectorAccessor<FFTWComplex> default_accessor; 00684 typedef VectorAccessor<FFTWComplex> DefaultAccessor; 00685 typedef VigraTrueType hasConstantStrides; 00686 }; 00687 00688 /** Complex (FFTWComplex) image. 00689 It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and 00690 their const counterparts to access the data. 00691 00692 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00693 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00694 Namespace: vigra 00695 */ 00696 typedef BasicImage<FFTWComplex> FFTWComplexImage; 00697 00698 //@} 00699 00700 /********************************************************/ 00701 /* */ 00702 /* FFTWComplex-Accessors */ 00703 /* */ 00704 /********************************************************/ 00705 00706 /** \addtogroup DataAccessors 00707 */ 00708 //@{ 00709 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex 00710 00711 Encapsulate access to pixels of type FFTWComplex 00712 */ 00713 //@{ 00714 /** Encapsulate access to the the real part of a complex number. 00715 00716 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00717 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00718 Namespace: vigra 00719 */ 00720 class FFTWRealAccessor 00721 { 00722 public: 00723 00724 /// The accessor's value type. 00725 typedef fftw_real value_type; 00726 00727 /// Read real part at iterator position. 00728 template <class ITERATOR> 00729 value_type operator()(ITERATOR const & i) const { 00730 return (*i).re(); 00731 } 00732 00733 /// Read real part at offset from iterator position. 00734 template <class ITERATOR, class DIFFERENCE> 00735 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00736 return i[d].re(); 00737 } 00738 00739 /// Write real part at iterator position from a scalar. 00740 template <class ITERATOR> 00741 void set(value_type const & v, ITERATOR const & i) const { 00742 (*i).re()= v; 00743 } 00744 00745 /// Write real part at offset from iterator position from a scalar. 00746 template <class ITERATOR, class DIFFERENCE> 00747 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 00748 i[d].re()= v; 00749 } 00750 00751 /// Write real part at iterator position into a scalar. 00752 template <class ITERATOR> 00753 void set(FFTWComplex const & v, ITERATOR const & i) const { 00754 *i = v.re(); 00755 } 00756 00757 /// Write real part at offset from iterator position into a scalar. 00758 template <class ITERATOR, class DIFFERENCE> 00759 void set(FFTWComplex const & v, ITERATOR const & i, DIFFERENCE d) const { 00760 i[d] = v.re(); 00761 } 00762 }; 00763 00764 /** Encapsulate access to the the imaginary part of a complex number. 00765 00766 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00767 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00768 Namespace: vigra 00769 */ 00770 class FFTWImaginaryAccessor 00771 { 00772 public: 00773 /// The accessor's value type. 00774 typedef fftw_real value_type; 00775 00776 /// Read imaginary part at iterator position. 00777 template <class ITERATOR> 00778 value_type operator()(ITERATOR const & i) const { 00779 return (*i).im(); 00780 } 00781 00782 /// Read imaginary part at offset from iterator position. 00783 template <class ITERATOR, class DIFFERENCE> 00784 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00785 return i[d].im(); 00786 } 00787 00788 /// Write imaginary part at iterator position from a scalar. 00789 template <class ITERATOR> 00790 void set(value_type const & v, ITERATOR const & i) const { 00791 (*i).im()= v; 00792 } 00793 00794 /// Write imaginary part at offset from iterator position from a scalar. 00795 template <class ITERATOR, class DIFFERENCE> 00796 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 00797 i[d].im()= v; 00798 } 00799 00800 /// Write imaginary part at iterator position into a scalar. 00801 template <class ITERATOR> 00802 void set(FFTWComplex const & v, ITERATOR const & i) const { 00803 *i = v.im(); 00804 } 00805 00806 /// Write imaginary part at offset from iterator position into a scalar. 00807 template <class ITERATOR, class DIFFERENCE> 00808 void set(FFTWComplex const & v, ITERATOR const & i, DIFFERENCE d) const { 00809 i[d] = v.im(); 00810 } 00811 }; 00812 00813 /** Write a real number into a complex one. The imaginary part is set 00814 to 0. 00815 00816 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00817 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00818 Namespace: vigra 00819 */ 00820 class FFTWWriteRealAccessor: public FFTWRealAccessor 00821 { 00822 public: 00823 /// The accessor's value type. 00824 typedef fftw_real value_type; 00825 00826 /** Write real number at iterator position. Set imaginary part 00827 to 0. 00828 */ 00829 template <class ITERATOR> 00830 void set(value_type const & v, ITERATOR const & i) const { 00831 (*i).re()= v; 00832 (*i).im()= 0; 00833 } 00834 00835 /** Write real number at offset from iterator position. Set imaginary part 00836 to 0. 00837 */ 00838 template <class ITERATOR, class DIFFERENCE> 00839 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 00840 i[d].re()= v; 00841 i[d].im()= 0; 00842 } 00843 }; 00844 00845 /** Calculate magnitude of complex number on the fly. 00846 00847 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00848 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00849 Namespace: vigra 00850 */ 00851 class FFTWMagnitudeAccessor 00852 { 00853 public: 00854 /// The accessor's value type. 00855 typedef fftw_real value_type; 00856 00857 /// Read magnitude at iterator position. 00858 template <class ITERATOR> 00859 value_type operator()(ITERATOR const & i) const { 00860 return (*i).magnitude(); 00861 } 00862 00863 /// Read magnitude at offset from iterator position. 00864 template <class ITERATOR, class DIFFERENCE> 00865 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00866 return (i[d]).magnitude(); 00867 } 00868 }; 00869 00870 /** Calculate phase of complex number on the fly. 00871 00872 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>> (for FFTW 3) or<br> 00873 <b>\#include</b> <<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>> (for deprecated FFTW 2)<br> 00874 Namespace: vigra 00875 */ 00876 class FFTWPhaseAccessor 00877 { 00878 public: 00879 /// The accessor's value type. 00880 typedef fftw_real value_type; 00881 00882 /// Read phase at iterator position. 00883 template <class ITERATOR> 00884 value_type operator()(ITERATOR const & i) const { 00885 return (*i).phase(); 00886 } 00887 00888 /// Read phase at offset from iterator position. 00889 template <class ITERATOR, class DIFFERENCE> 00890 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00891 return (i[d]).phase(); 00892 } 00893 }; 00894 00895 //@} 00896 //@} 00897 00898 /********************************************************/ 00899 /* */ 00900 /* Fourier Transform */ 00901 /* */ 00902 /********************************************************/ 00903 00904 /** \addtogroup FourierTransform Fast Fourier Transform 00905 00906 This documentation describes the VIGRA interface to FFTW version 3. The interface 00907 to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated. 00908 00909 VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier 00910 Transform</a> package to perform Fourier transformations. VIGRA 00911 provides a wrapper for FFTW's complex number type (FFTWComplex), 00912 but FFTW's functions are used verbatim. If the image is stored as 00913 a FFTWComplexImage, the simplest call to an FFT function is like this: 00914 00915 \code 00916 vigra::FFTWComplexImage spatial(width,height), fourier(width,height); 00917 ... // fill image with data 00918 00919 // create a plan with estimated performance optimization 00920 fftw_plan forwardPlan = fftw_plan_dft_2d(height, width, 00921 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), 00922 FFTW_FORWARD, FFTW_ESTIMATE ); 00923 // calculate FFT (this can be repeated as often as needed, 00924 // with fresh data written into the source array) 00925 fftw_execute(forwardPlan); 00926 00927 // release the plan memory 00928 fftw_destroy_plan(forwardPlan); 00929 00930 // likewise for the inverse transform 00931 fftw_plan backwardPlan = fftw_plan_dft_2d(height, width, 00932 (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(), 00933 FFTW_BACKWARD, FFTW_ESTIMATE); 00934 fftw_execute(backwardPlan); 00935 fftw_destroy_plan(backwardPlan); 00936 00937 // do not forget to normalize the result according to the image size 00938 transformImage(srcImageRange(spatial), destImage(spatial), 00939 std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height)); 00940 \endcode 00941 00942 Note that in the creation of a plan, the height must be given 00943 first. Note also that <TT>spatial.begin()</TT> may only be passed 00944 to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the 00945 entire image. When you want to restrict operation to an ROI, you 00946 can create a copy of the ROI in an image of appropriate size, or 00947 you may use the Guru interface to FFTW. 00948 00949 More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>. 00950 00951 FFTW produces fourier images that have the DC component (the 00952 origin of the Fourier space) in the upper left corner. Often, one 00953 wants the origin in the center of the image, so that frequencies 00954 always increase towards the border of the image. This can be 00955 achieved by calling \ref moveDCToCenter(). The inverse 00956 transformation is done by \ref moveDCToUpperLeft(). 00957 00958 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br> 00959 Namespace: vigra 00960 */ 00961 00962 /** \addtogroup FourierTransform 00963 */ 00964 //@{ 00965 00966 /********************************************************/ 00967 /* */ 00968 /* moveDCToCenter */ 00969 /* */ 00970 /********************************************************/ 00971 00972 /** \brief Rearrange the quadrants of a Fourier image so that the origin is 00973 in the image center. 00974 00975 FFTW produces fourier images where the DC component (origin of 00976 fourier space) is located in the upper left corner of the 00977 image. The quadrants are placed like this (using a 4x4 image for 00978 example): 00979 00980 \code 00981 DC 4 3 3 00982 4 4 3 3 00983 1 1 2 2 00984 1 1 2 2 00985 \endcode 00986 00987 After applying the function, the quadrants are at their usual places: 00988 00989 \code 00990 2 2 1 1 00991 2 2 1 1 00992 3 3 DC 4 00993 3 3 4 4 00994 \endcode 00995 00996 This transformation can be reversed by \ref moveDCToUpperLeft(). 00997 Note that the transformation must not be executed in place - input 00998 and output images must be different. 00999 01000 <b> Declarations:</b> 01001 01002 pass arguments explicitly: 01003 \code 01004 namespace vigra { 01005 template <class SrcImageIterator, class SrcAccessor, 01006 class DestImageIterator, class DestAccessor> 01007 void moveDCToCenter(SrcImageIterator src_upperleft, 01008 SrcImageIterator src_lowerright, SrcAccessor sa, 01009 DestImageIterator dest_upperleft, DestAccessor da); 01010 } 01011 \endcode 01012 01013 01014 use argument objects in conjunction with \ref ArgumentObjectFactories : 01015 \code 01016 namespace vigra { 01017 template <class SrcImageIterator, class SrcAccessor, 01018 class DestImageIterator, class DestAccessor> 01019 void moveDCToCenter( 01020 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01021 pair<DestImageIterator, DestAccessor> dest); 01022 } 01023 \endcode 01024 01025 <b> Usage:</b> 01026 01027 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br> 01028 Namespace: vigra 01029 01030 \code 01031 vigra::FFTWComplexImage spatial(width,height), fourier(width,height); 01032 ... // fill image with data 01033 01034 // create a plan with estimated performance optimization 01035 fftw_plan forwardPlan = fftw_plan_dft_2d(height, width, 01036 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), 01037 FFTW_FORWARD, FFTW_ESTIMATE ); 01038 // calculate FFT 01039 fftw_execute(forwardPlan); 01040 01041 vigra::FFTWComplexImage rearrangedFourier(width, height); 01042 moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier)); 01043 01044 // delete the plan 01045 fftw_destroy_plan(forwardPlan); 01046 \endcode 01047 */ 01048 doxygen_overloaded_function(template <...> void moveDCToCenter) 01049 01050 template <class SrcImageIterator, class SrcAccessor, 01051 class DestImageIterator, class DestAccessor> 01052 void moveDCToCenter(SrcImageIterator src_upperleft, 01053 SrcImageIterator src_lowerright, SrcAccessor sa, 01054 DestImageIterator dest_upperleft, DestAccessor da) 01055 { 01056 int w = int(src_lowerright.x - src_upperleft.x); 01057 int h = int(src_lowerright.y - src_upperleft.y); 01058 int w1 = w/2; 01059 int h1 = h/2; 01060 int w2 = (w+1)/2; 01061 int h2 = (h+1)/2; 01062 01063 // 2. Quadrant zum 4. 01064 copyImage(srcIterRange(src_upperleft, 01065 src_upperleft + Diff2D(w2, h2), sa), 01066 destIter (dest_upperleft + Diff2D(w1, h1), da)); 01067 01068 // 4. Quadrant zum 2. 01069 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2), 01070 src_lowerright, sa), 01071 destIter (dest_upperleft, da)); 01072 01073 // 1. Quadrant zum 3. 01074 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0), 01075 src_upperleft + Diff2D(w, h2), sa), 01076 destIter (dest_upperleft + Diff2D(0, h1), da)); 01077 01078 // 3. Quadrant zum 1. 01079 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2), 01080 src_upperleft + Diff2D(w2, h), sa), 01081 destIter (dest_upperleft + Diff2D(w1, 0), da)); 01082 } 01083 01084 template <class SrcImageIterator, class SrcAccessor, 01085 class DestImageIterator, class DestAccessor> 01086 inline void moveDCToCenter( 01087 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01088 pair<DestImageIterator, DestAccessor> dest) 01089 { 01090 moveDCToCenter(src.first, src.second, src.third, 01091 dest.first, dest.second); 01092 } 01093 01094 /********************************************************/ 01095 /* */ 01096 /* moveDCToUpperLeft */ 01097 /* */ 01098 /********************************************************/ 01099 01100 /** \brief Rearrange the quadrants of a Fourier image so that the origin is 01101 in the image's upper left. 01102 01103 This function is the inversion of \ref moveDCToCenter(). See there 01104 for declarations and a usage example. 01105 01106 <b> Declarations:</b> 01107 01108 pass arguments explicitly: 01109 \code 01110 namespace vigra { 01111 template <class SrcImageIterator, class SrcAccessor, 01112 class DestImageIterator, class DestAccessor> 01113 void moveDCToUpperLeft(SrcImageIterator src_upperleft, 01114 SrcImageIterator src_lowerright, SrcAccessor sa, 01115 DestImageIterator dest_upperleft, DestAccessor da); 01116 } 01117 \endcode 01118 01119 01120 use argument objects in conjunction with \ref ArgumentObjectFactories : 01121 \code 01122 namespace vigra { 01123 template <class SrcImageIterator, class SrcAccessor, 01124 class DestImageIterator, class DestAccessor> 01125 void moveDCToUpperLeft( 01126 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01127 pair<DestImageIterator, DestAccessor> dest); 01128 } 01129 \endcode 01130 */ 01131 doxygen_overloaded_function(template <...> void moveDCToUpperLeft) 01132 01133 template <class SrcImageIterator, class SrcAccessor, 01134 class DestImageIterator, class DestAccessor> 01135 void moveDCToUpperLeft(SrcImageIterator src_upperleft, 01136 SrcImageIterator src_lowerright, SrcAccessor sa, 01137 DestImageIterator dest_upperleft, DestAccessor da) 01138 { 01139 int w = int(src_lowerright.x - src_upperleft.x); 01140 int h = int(src_lowerright.y - src_upperleft.y); 01141 int w2 = w/2; 01142 int h2 = h/2; 01143 int w1 = (w+1)/2; 01144 int h1 = (h+1)/2; 01145 01146 // 2. Quadrant zum 4. 01147 copyImage(srcIterRange(src_upperleft, 01148 src_upperleft + Diff2D(w2, h2), sa), 01149 destIter (dest_upperleft + Diff2D(w1, h1), da)); 01150 01151 // 4. Quadrant zum 2. 01152 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2), 01153 src_lowerright, sa), 01154 destIter (dest_upperleft, da)); 01155 01156 // 1. Quadrant zum 3. 01157 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0), 01158 src_upperleft + Diff2D(w, h2), sa), 01159 destIter (dest_upperleft + Diff2D(0, h1), da)); 01160 01161 // 3. Quadrant zum 1. 01162 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2), 01163 src_upperleft + Diff2D(w2, h), sa), 01164 destIter (dest_upperleft + Diff2D(w1, 0), da)); 01165 } 01166 01167 template <class SrcImageIterator, class SrcAccessor, 01168 class DestImageIterator, class DestAccessor> 01169 inline void moveDCToUpperLeft( 01170 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01171 pair<DestImageIterator, DestAccessor> dest) 01172 { 01173 moveDCToUpperLeft(src.first, src.second, src.third, 01174 dest.first, dest.second); 01175 } 01176 01177 template <class DestImageIterator, class DestAccessor> 01178 void fftShift(DestImageIterator upperleft, 01179 DestImageIterator lowerright, DestAccessor da) 01180 { 01181 int w = int(lowerright.x - upperleft.x); 01182 int h = int(lowerright.y - upperleft.y); 01183 int w2 = w/2; 01184 int h2 = h/2; 01185 int w1 = (w+1)/2; 01186 int h1 = (h+1)/2; 01187 01188 // 2. Quadrant zum 4. 01189 swapImageData(destIterRange(upperleft, 01190 upperleft + Diff2D(w2, h2), da), 01191 destIter (upperleft + Diff2D(w1, h1), da)); 01192 01193 // 1. Quadrant zum 3. 01194 swapImageData(destIterRange(upperleft + Diff2D(w2, 0), 01195 upperleft + Diff2D(w, h2), da), 01196 destIter (upperleft + Diff2D(0, h1), da)); 01197 } 01198 01199 template <class DestImageIterator, class DestAccessor> 01200 inline void fftShift( 01201 triple<DestImageIterator, DestImageIterator, DestAccessor> dest) 01202 { 01203 fftShift(dest.first, dest.second, dest.third); 01204 } 01205 01206 01207 namespace detail { 01208 01209 template <class T> 01210 void 01211 fourierTransformImpl(FFTWComplexImage::const_traverser sul, 01212 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01213 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest, T sign) 01214 { 01215 int w = int(slr.x - sul.x); 01216 int h = int(slr.y - sul.y); 01217 01218 FFTWComplexImage sworkImage, dworkImage; 01219 01220 fftw_complex * srcPtr = (fftw_complex *)(&*sul); 01221 fftw_complex * destPtr = (fftw_complex *)(&*dul); 01222 01223 // test for right memory layout (fftw expects a 2*width*height floats array) 01224 if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1)))) 01225 { 01226 sworkImage.resize(w, h); 01227 copyImage(srcIterRange(sul, slr, src), destImage(sworkImage)); 01228 srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft())); 01229 } 01230 if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1)))) 01231 { 01232 dworkImage.resize(w, h); 01233 destPtr = (fftw_complex *)(&(*dworkImage.upperLeft())); 01234 } 01235 01236 fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE ); 01237 fftw_execute(plan); 01238 fftw_destroy_plan(plan); 01239 01240 if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1)))) 01241 { 01242 copyImage(srcImageRange(dworkImage), destIter(dul, dest)); 01243 } 01244 } 01245 01246 } // namespace detail 01247 01248 /********************************************************/ 01249 /* */ 01250 /* fourierTransform */ 01251 /* */ 01252 /********************************************************/ 01253 01254 /** \brief Compute forward and inverse Fourier transforms. 01255 01256 In the forward direction, the input image may be scalar or complex, and the output image 01257 is always complex. In the inverse direction, both input and output must be complex. 01258 01259 <b> Declarations:</b> 01260 01261 pass arguments explicitly: 01262 \code 01263 namespace vigra { 01264 template <class SrcImageIterator, class SrcAccessor> 01265 void fourierTransform(SrcImageIterator srcUpperLeft, 01266 SrcImageIterator srcLowerRight, SrcAccessor src, 01267 FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest); 01268 01269 void 01270 fourierTransformInverse(FFTWComplexImage::const_traverser sul, 01271 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01272 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest) 01273 } 01274 \endcode 01275 01276 use argument objects in conjunction with \ref ArgumentObjectFactories : 01277 \code 01278 namespace vigra { 01279 template <class SrcImageIterator, class SrcAccessor> 01280 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01281 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest); 01282 01283 void 01284 fourierTransformInverse(triple<FFTWComplexImage::const_traverser, 01285 FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src, 01286 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest); 01287 } 01288 \endcode 01289 01290 <b> Usage:</b> 01291 01292 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br> 01293 Namespace: vigra 01294 01295 \code 01296 // compute complex Fourier transform of a real image 01297 vigra::DImage src(w, h); 01298 vigra::FFTWComplexImage fourier(w, h); 01299 01300 fourierTransform(srcImageRange(src), destImage(fourier)); 01301 01302 // compute inverse Fourier transform 01303 // note that both source and destination image must be of type vigra::FFTWComplexImage 01304 vigra::FFTWComplexImage inverseFourier(w, h); 01305 01306 fourierTransform(srcImageRange(fourier), destImage(inverseFourier)); 01307 \endcode 01308 */ 01309 doxygen_overloaded_function(template <...> void fourierTransform) 01310 01311 inline void 01312 fourierTransform(FFTWComplexImage::const_traverser sul, 01313 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01314 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest) 01315 { 01316 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD); 01317 } 01318 01319 template <class SrcImageIterator, class SrcAccessor> 01320 void fourierTransform(SrcImageIterator srcUpperLeft, 01321 SrcImageIterator srcLowerRight, SrcAccessor sa, 01322 FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor da) 01323 { 01324 // copy real input images into a complex one... 01325 int w= srcLowerRight.x - srcUpperLeft.x; 01326 int h= srcLowerRight.y - srcUpperLeft.y; 01327 01328 FFTWComplexImage workImage(w, h); 01329 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01330 destImage(workImage, FFTWWriteRealAccessor())); 01331 01332 // ...and call the complex -> complex version of the algorithm 01333 FFTWComplexImage const & cworkImage = workImage; 01334 fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01335 destUpperLeft, da); 01336 } 01337 01338 template <class SrcImageIterator, class SrcAccessor> 01339 inline 01340 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01341 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest) 01342 { 01343 fourierTransform(src.first, src.second, src.third, dest.first, dest.second); 01344 } 01345 01346 /** \brief Compute inverse Fourier transforms. 01347 01348 See \ref fourierTransform() for details. 01349 */ 01350 inline void 01351 fourierTransformInverse(FFTWComplexImage::const_traverser sul, 01352 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01353 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest) 01354 { 01355 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD); 01356 } 01357 01358 template <class DestImageIterator, class DestAccessor> 01359 void fourierTransformInverse(FFTWComplexImage::const_traverser sul, 01360 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01361 DestImageIterator dul, DestAccessor dest) 01362 { 01363 int w = slr.x - sul.x; 01364 int h = slr.y - sul.y; 01365 01366 FFTWComplexImage workImage(w, h); 01367 fourierTransformInverse(sul, slr, src, workImage.upperLeft(), workImage.accessor()); 01368 copyImage(srcImageRange(workImage), destIter(dul, dest)); 01369 } 01370 01371 01372 template <class DestImageIterator, class DestAccessor> 01373 inline void 01374 fourierTransformInverse(triple<FFTWComplexImage::const_traverser, 01375 FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src, 01376 pair<DestImageIterator, DestAccessor> dest) 01377 { 01378 fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second); 01379 } 01380 01381 01382 01383 /********************************************************/ 01384 /* */ 01385 /* applyFourierFilter */ 01386 /* */ 01387 /********************************************************/ 01388 01389 /** \brief Apply a filter (defined in the frequency domain) to an image. 01390 01391 After transferring the image into the frequency domain, it is 01392 multiplied pixel-wise with the filter and transformed back. The 01393 result is put into the given destination image which must have the right size. 01394 The result will be normalized to compensate for the two FFTs. 01395 01396 If the destination image is scalar, only the real part of the result image is 01397 retained. In this case, you are responsible for choosing a filter image 01398 which ensures a zero imaginary part of the result (e.g. use a real, even symmetric 01399 filter image, or a purely imaginary, odd symmetric on). 01400 01401 The DC entry of the filter must be in the upper left, which is the 01402 position where FFTW expects it (see \ref moveDCToUpperLeft()). 01403 01404 <b> Declarations:</b> 01405 01406 pass arguments explicitly: 01407 \code 01408 namespace vigra { 01409 template <class SrcImageIterator, class SrcAccessor, 01410 class FilterImageIterator, class FilterAccessor, 01411 class DestImageIterator, class DestAccessor> 01412 void applyFourierFilter(SrcImageIterator srcUpperLeft, 01413 SrcImageIterator srcLowerRight, SrcAccessor sa, 01414 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01415 DestImageIterator destUpperLeft, DestAccessor da); 01416 } 01417 \endcode 01418 01419 use argument objects in conjunction with \ref ArgumentObjectFactories : 01420 \code 01421 namespace vigra { 01422 template <class SrcImageIterator, class SrcAccessor, 01423 class FilterImageIterator, class FilterAccessor, 01424 class DestImageIterator, class DestAccessor> 01425 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01426 pair<FilterImageIterator, FilterAccessor> filter, 01427 pair<DestImageIterator, DestAccessor> dest); 01428 } 01429 \endcode 01430 01431 <b> Usage:</b> 01432 01433 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br> 01434 Namespace: vigra 01435 01436 \code 01437 // create a Gaussian filter in Fourier space 01438 vigra::FImage gaussFilter(w, h), filter(w, h); 01439 for(int y=0; y<h; ++y) 01440 for(int x=0; x<w; ++x) 01441 { 01442 xx = float(x - w / 2) / w; 01443 yy = float(y - h / 2) / h; 01444 01445 gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale); 01446 } 01447 01448 // applyFourierFilter() expects the filter's DC in the upper left 01449 moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter)); 01450 01451 vigra::FFTWComplexImage result(w, h); 01452 01453 vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result); 01454 \endcode 01455 01456 For inspection of the result, \ref FFTWMagnitudeAccessor might be 01457 useful. If you want to apply the same filter repeatedly, it may be more 01458 efficient to use the FFTW functions directly with FFTW plans optimized 01459 for good performance. 01460 */ 01461 doxygen_overloaded_function(template <...> void applyFourierFilter) 01462 01463 template <class SrcImageIterator, class SrcAccessor, 01464 class FilterImageIterator, class FilterAccessor, 01465 class DestImageIterator, class DestAccessor> 01466 void applyFourierFilter(SrcImageIterator srcUpperLeft, 01467 SrcImageIterator srcLowerRight, SrcAccessor sa, 01468 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01469 DestImageIterator destUpperLeft, DestAccessor da) 01470 { 01471 // copy real input images into a complex one... 01472 int w = int(srcLowerRight.x - srcUpperLeft.x); 01473 int h = int(srcLowerRight.y - srcUpperLeft.y); 01474 01475 FFTWComplexImage workImage(w, h); 01476 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01477 destImage(workImage, FFTWWriteRealAccessor())); 01478 01479 // ...and call the impl 01480 FFTWComplexImage const & cworkImage = workImage; 01481 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01482 filterUpperLeft, fa, 01483 destUpperLeft, da); 01484 } 01485 01486 template <class FilterImageIterator, class FilterAccessor, 01487 class DestImageIterator, class DestAccessor> 01488 inline 01489 void applyFourierFilter( 01490 FFTWComplexImage::const_traverser srcUpperLeft, 01491 FFTWComplexImage::const_traverser srcLowerRight, 01492 FFTWComplexImage::ConstAccessor sa, 01493 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01494 DestImageIterator destUpperLeft, DestAccessor da) 01495 { 01496 int w = srcLowerRight.x - srcUpperLeft.x; 01497 int h = srcLowerRight.y - srcUpperLeft.y; 01498 01499 // test for right memory layout (fftw expects a 2*width*height floats array) 01500 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 01501 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa, 01502 filterUpperLeft, fa, 01503 destUpperLeft, da); 01504 else 01505 { 01506 FFTWComplexImage workImage(w, h); 01507 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01508 destImage(workImage)); 01509 01510 FFTWComplexImage const & cworkImage = workImage; 01511 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01512 filterUpperLeft, fa, 01513 destUpperLeft, da); 01514 } 01515 } 01516 01517 template <class SrcImageIterator, class SrcAccessor, 01518 class FilterImageIterator, class FilterAccessor, 01519 class DestImageIterator, class DestAccessor> 01520 inline 01521 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01522 pair<FilterImageIterator, FilterAccessor> filter, 01523 pair<DestImageIterator, DestAccessor> dest) 01524 { 01525 applyFourierFilter(src.first, src.second, src.third, 01526 filter.first, filter.second, 01527 dest.first, dest.second); 01528 } 01529 01530 template <class FilterImageIterator, class FilterAccessor, 01531 class DestImageIterator, class DestAccessor> 01532 void applyFourierFilterImpl( 01533 FFTWComplexImage::const_traverser srcUpperLeft, 01534 FFTWComplexImage::const_traverser srcLowerRight, 01535 FFTWComplexImage::ConstAccessor, 01536 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01537 DestImageIterator destUpperLeft, DestAccessor da) 01538 { 01539 int w = int(srcLowerRight.x - srcUpperLeft.x); 01540 int h = int(srcLowerRight.y - srcUpperLeft.y); 01541 01542 FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft); 01543 01544 // FFT from srcImage to complexResultImg 01545 fftw_plan forwardPlan= 01546 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft), 01547 (fftw_complex *)complexResultImg.begin(), 01548 FFTW_FORWARD, FFTW_ESTIMATE ); 01549 fftw_execute(forwardPlan); 01550 fftw_destroy_plan(forwardPlan); 01551 01552 // convolve in freq. domain (in complexResultImg) 01553 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa), 01554 destImage(complexResultImg), std::multiplies<FFTWComplex>()); 01555 01556 // FFT back into spatial domain (inplace in complexResultImg) 01557 fftw_plan backwardPlan= 01558 fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(), 01559 (fftw_complex *)complexResultImg.begin(), 01560 FFTW_BACKWARD, FFTW_ESTIMATE); 01561 fftw_execute(backwardPlan); 01562 fftw_destroy_plan(backwardPlan); 01563 01564 typedef typename 01565 NumericTraits<typename DestAccessor::value_type>::isScalar 01566 isScalarResult; 01567 01568 // normalization (after FFTs), maybe stripping imaginary part 01569 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da, 01570 isScalarResult()); 01571 } 01572 01573 template <class DestImageIterator, class DestAccessor> 01574 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage, 01575 DestImageIterator destUpperLeft, 01576 DestAccessor da, 01577 VigraFalseType) 01578 { 01579 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 01580 01581 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 01582 { 01583 DestImageIterator dIt= destUpperLeft; 01584 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 01585 { 01586 da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0); 01587 da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1); 01588 } 01589 } 01590 } 01591 01592 inline 01593 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 01594 FFTWComplexImage::traverser destUpperLeft, 01595 FFTWComplexImage::Accessor da, 01596 VigraFalseType) 01597 { 01598 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da), 01599 linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height()))); 01600 } 01601 01602 template <class DestImageIterator, class DestAccessor> 01603 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 01604 DestImageIterator destUpperLeft, 01605 DestAccessor da, 01606 VigraTrueType) 01607 { 01608 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 01609 01610 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 01611 { 01612 DestImageIterator dIt= destUpperLeft; 01613 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 01614 da.set(srcImage(x, y).re()*normFactor, dIt); 01615 } 01616 } 01617 01618 /**********************************************************/ 01619 /* */ 01620 /* applyFourierFilterFamily */ 01621 /* */ 01622 /**********************************************************/ 01623 01624 /** \brief Apply an array of filters (defined in the frequency domain) to an image. 01625 01626 This provides the same functionality as \ref applyFourierFilter(), 01627 but applying several filters at once allows to avoid 01628 repeated Fourier transforms of the source image. 01629 01630 Filters and result images must be stored in \ref vigra::ImageArray data 01631 structures. In contrast to \ref applyFourierFilter(), this function adjusts 01632 the size of the result images and the the length of the array. 01633 01634 <b> Declarations:</b> 01635 01636 pass arguments explicitly: 01637 \code 01638 namespace vigra { 01639 template <class SrcImageIterator, class SrcAccessor, class FilterType> 01640 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 01641 SrcImageIterator srcLowerRight, SrcAccessor sa, 01642 const ImageArray<FilterType> &filters, 01643 ImageArray<FFTWComplexImage> &results) 01644 } 01645 \endcode 01646 01647 use argument objects in conjunction with \ref ArgumentObjectFactories : 01648 \code 01649 namespace vigra { 01650 template <class SrcImageIterator, class SrcAccessor, class FilterType> 01651 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01652 const ImageArray<FilterType> &filters, 01653 ImageArray<FFTWComplexImage> &results) 01654 } 01655 \endcode 01656 01657 <b> Usage:</b> 01658 01659 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br> 01660 Namespace: vigra 01661 01662 \code 01663 // assuming the presence of a real-valued image named "image" to 01664 // be filtered in this example 01665 01666 vigra::ImageArray<vigra::FImage> filters(16, image.size()); 01667 01668 for (int i=0; i<filters.size(); i++) 01669 // create some meaningful filters here 01670 createMyFilterOfScale(i, destImage(filters[i])); 01671 01672 vigra::ImageArray<vigra::FFTWComplexImage> results(); 01673 01674 vigra::applyFourierFilterFamily(srcImageRange(image), filters, results); 01675 \endcode 01676 */ 01677 doxygen_overloaded_function(template <...> void applyFourierFilterFamily) 01678 01679 template <class SrcImageIterator, class SrcAccessor, 01680 class FilterType, class DestImage> 01681 inline 01682 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01683 const ImageArray<FilterType> &filters, 01684 ImageArray<DestImage> &results) 01685 { 01686 applyFourierFilterFamily(src.first, src.second, src.third, 01687 filters, results); 01688 } 01689 01690 template <class SrcImageIterator, class SrcAccessor, 01691 class FilterType, class DestImage> 01692 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 01693 SrcImageIterator srcLowerRight, SrcAccessor sa, 01694 const ImageArray<FilterType> &filters, 01695 ImageArray<DestImage> &results) 01696 { 01697 int w = int(srcLowerRight.x - srcUpperLeft.x); 01698 int h = int(srcLowerRight.y - srcUpperLeft.y); 01699 01700 FFTWComplexImage workImage(w, h); 01701 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01702 destImage(workImage, FFTWWriteRealAccessor())); 01703 01704 FFTWComplexImage const & cworkImage = workImage; 01705 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01706 filters, results); 01707 } 01708 01709 template <class FilterType, class DestImage> 01710 inline 01711 void applyFourierFilterFamily( 01712 FFTWComplexImage::const_traverser srcUpperLeft, 01713 FFTWComplexImage::const_traverser srcLowerRight, 01714 FFTWComplexImage::ConstAccessor sa, 01715 const ImageArray<FilterType> &filters, 01716 ImageArray<DestImage> &results) 01717 { 01718 int w= srcLowerRight.x - srcUpperLeft.x; 01719 01720 // test for right memory layout (fftw expects a 2*width*height floats array) 01721 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 01722 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa, 01723 filters, results); 01724 else 01725 { 01726 int h = srcLowerRight.y - srcUpperLeft.y; 01727 FFTWComplexImage workImage(w, h); 01728 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01729 destImage(workImage)); 01730 01731 FFTWComplexImage const & cworkImage = workImage; 01732 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01733 filters, results); 01734 } 01735 } 01736 01737 template <class FilterType, class DestImage> 01738 void applyFourierFilterFamilyImpl( 01739 FFTWComplexImage::const_traverser srcUpperLeft, 01740 FFTWComplexImage::const_traverser srcLowerRight, 01741 FFTWComplexImage::ConstAccessor sa, 01742 const ImageArray<FilterType> &filters, 01743 ImageArray<DestImage> &results) 01744 { 01745 // FIXME: sa is not used 01746 // (maybe check if StandardAccessor, else copy?) 01747 01748 // make sure the filter images have the right dimensions 01749 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(), 01750 "applyFourierFilterFamily called with src image size != filters.imageSize()!"); 01751 01752 // make sure the result image array has the right dimensions 01753 results.resize(filters.size()); 01754 results.resizeImages(filters.imageSize()); 01755 01756 // FFT from srcImage to freqImage 01757 int w = int(srcLowerRight.x - srcUpperLeft.x); 01758 int h = int(srcLowerRight.y - srcUpperLeft.y); 01759 01760 FFTWComplexImage freqImage(w, h); 01761 FFTWComplexImage result(w, h); 01762 01763 fftw_plan forwardPlan= 01764 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft), 01765 (fftw_complex *)freqImage.begin(), 01766 FFTW_FORWARD, FFTW_ESTIMATE ); 01767 fftw_execute(forwardPlan); 01768 fftw_destroy_plan(forwardPlan); 01769 01770 fftw_plan backwardPlan= 01771 fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(), 01772 (fftw_complex *)result.begin(), 01773 FFTW_BACKWARD, FFTW_ESTIMATE ); 01774 typedef typename 01775 NumericTraits<typename DestImage::Accessor::value_type>::isScalar 01776 isScalarResult; 01777 01778 // convolve with filters in freq. domain 01779 for (unsigned int i= 0; i < filters.size(); i++) 01780 { 01781 combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]), 01782 destImage(result), std::multiplies<FFTWComplex>()); 01783 01784 // FFT back into spatial domain (inplace in destImage) 01785 fftw_execute(backwardPlan); 01786 01787 // normalization (after FFTs), maybe stripping imaginary part 01788 applyFourierFilterImplNormalization(result, 01789 results[i].upperLeft(), results[i].accessor(), 01790 isScalarResult()); 01791 } 01792 fftw_destroy_plan(backwardPlan); 01793 } 01794 01795 /********************************************************/ 01796 /* */ 01797 /* fourierTransformReal */ 01798 /* */ 01799 /********************************************************/ 01800 01801 /** \brief Real Fourier transforms for even and odd boundary conditions 01802 (aka. cosine and sine transforms). 01803 01804 01805 If the image is real and has even symmetry, its Fourier transform 01806 is also real and has even symmetry. The Fourier transform of a real image with odd 01807 symmetry is imaginary and has odd symmetry. In either case, only about a quarter 01808 of the pixels need to be stored because the rest can be calculated from the symmetry 01809 properties. This is especially useful, if the original image is implicitly assumed 01810 to have reflective or anti-reflective boundary conditions. Then the "negative" 01811 pixel locations are defined as 01812 01813 \code 01814 even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1) 01815 odd (anti-reflective boundary conditions): f[-1] = 0 01816 f[-x] = -f[x-2] (x = 2,...,N-1) 01817 \endcode 01818 01819 end similar at the other boundary (see the FFTW documentation for details). 01820 This has the advantage that more efficient Fourier transforms that use only 01821 real numbers can be implemented. These are also known as cosine and sine transforms 01822 respectively. 01823 01824 If you use the odd transform it is important to note that in the Fourier domain, 01825 the DC component is always zero and is therefore dropped from the data structure. 01826 This means that index 0 in an odd symmetric Fourier domain image refers to 01827 the <i>first</i> harmonic. This is especially important if an image is first 01828 cosine transformed (even symmetry), then in the Fourier domain multiplied 01829 with an odd symmetric filter (e.g. a first derivative) and finally transformed 01830 back to the spatial domain with a sine transform (odd symmetric). For this to work 01831 properly the image must be shifted left or up by one pixel (depending on whether 01832 the x- or y-axis is odd symmetric) before the inverse transform can be applied. 01833 (see example below). 01834 01835 The real Fourier transform functions are named <tt>fourierTransformReal??</tt> 01836 where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating 01837 whether the x- and y-axis is to be transformed using even or odd symmetry. 01838 The same functions can be used for both the forward and inverse transforms, 01839 only the normalization changes. For signal processing, the following 01840 normalization factors are most appropriate: 01841 01842 \code 01843 forward inverse 01844 ------------------------------------------------------------ 01845 X even, Y even 1.0 4.0 * (w-1) * (h-1) 01846 X even, Y odd -1.0 -4.0 * (w-1) * (h+1) 01847 X odd, Y even -1.0 -4.0 * (w+1) * (h-1) 01848 X odd, Y odd 1.0 4.0 * (w+1) * (h+1) 01849 \endcode 01850 01851 where <tt>w</tt> and <tt>h</tt> denote the image width and height. 01852 01853 <b> Declarations:</b> 01854 01855 pass arguments explicitly: 01856 \code 01857 namespace vigra { 01858 template <class SrcTraverser, class SrcAccessor, 01859 class DestTraverser, class DestAccessor> 01860 void 01861 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01862 DestTraverser dul, DestAccessor dest, fftw_real norm); 01863 01864 fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise 01865 } 01866 \endcode 01867 01868 01869 use argument objects in conjunction with \ref ArgumentObjectFactories : 01870 \code 01871 namespace vigra { 01872 template <class SrcTraverser, class SrcAccessor, 01873 class DestTraverser, class DestAccessor> 01874 void 01875 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01876 pair<DestTraverser, DestAccessor> dest, fftw_real norm); 01877 01878 fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise 01879 } 01880 \endcode 01881 01882 <b> Usage:</b> 01883 01884 <b>\#include</b> <<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>><br> 01885 Namespace: vigra 01886 01887 \code 01888 vigra::FImage spatial(width,height), fourier(width,height); 01889 ... // fill image with data 01890 01891 // forward cosine transform == reflective boundary conditions 01892 fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0); 01893 01894 // multiply with a first derivative of Gaussian in x-direction 01895 for(int y = 0; y < height; ++y) 01896 { 01897 for(int x = 1; x < width; ++x) 01898 { 01899 double dx = x * M_PI / (width - 1); 01900 double dy = y * M_PI / (height - 1); 01901 fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0); 01902 } 01903 fourier(width-1, y) = 0.0; 01904 } 01905 01906 // inverse transform -- odd symmetry in x-direction, even in y, 01907 // due to symmetry of the filter 01908 fourierTransformRealOE(srcImageRange(fourier), destImage(spatial), 01909 (fftw_real)-4.0 * (width+1) * (height-1)); 01910 \endcode 01911 */ 01912 doxygen_overloaded_function(template <...> void fourierTransformReal) 01913 01914 template <class SrcTraverser, class SrcAccessor, 01915 class DestTraverser, class DestAccessor> 01916 inline void 01917 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01918 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 01919 { 01920 fourierTransformRealEE(src.first, src.second, src.third, 01921 dest.first, dest.second, norm); 01922 } 01923 01924 template <class SrcTraverser, class SrcAccessor, 01925 class DestTraverser, class DestAccessor> 01926 inline void 01927 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01928 DestTraverser dul, DestAccessor dest, fftw_real norm) 01929 { 01930 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01931 norm, FFTW_REDFT00, FFTW_REDFT00); 01932 } 01933 01934 template <class DestTraverser, class DestAccessor> 01935 inline void 01936 fourierTransformRealEE( 01937 FFTWRealImage::const_traverser sul, 01938 FFTWRealImage::const_traverser slr, 01939 FFTWRealImage::Accessor src, 01940 DestTraverser dul, DestAccessor dest, fftw_real norm) 01941 { 01942 int w = slr.x - sul.x; 01943 01944 // test for right memory layout (fftw expects a width*height fftw_real array) 01945 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 01946 fourierTransformRealImpl(sul, slr, dul, dest, 01947 norm, FFTW_REDFT00, FFTW_REDFT00); 01948 else 01949 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01950 norm, FFTW_REDFT00, FFTW_REDFT00); 01951 } 01952 01953 /********************************************************************/ 01954 01955 template <class SrcTraverser, class SrcAccessor, 01956 class DestTraverser, class DestAccessor> 01957 inline void 01958 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01959 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 01960 { 01961 fourierTransformRealOE(src.first, src.second, src.third, 01962 dest.first, dest.second, norm); 01963 } 01964 01965 template <class SrcTraverser, class SrcAccessor, 01966 class DestTraverser, class DestAccessor> 01967 inline void 01968 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01969 DestTraverser dul, DestAccessor dest, fftw_real norm) 01970 { 01971 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01972 norm, FFTW_RODFT00, FFTW_REDFT00); 01973 } 01974 01975 template <class DestTraverser, class DestAccessor> 01976 inline void 01977 fourierTransformRealOE( 01978 FFTWRealImage::const_traverser sul, 01979 FFTWRealImage::const_traverser slr, 01980 FFTWRealImage::Accessor src, 01981 DestTraverser dul, DestAccessor dest, fftw_real norm) 01982 { 01983 int w = slr.x - sul.x; 01984 01985 // test for right memory layout (fftw expects a width*height fftw_real array) 01986 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 01987 fourierTransformRealImpl(sul, slr, dul, dest, 01988 norm, FFTW_RODFT00, FFTW_REDFT00); 01989 else 01990 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01991 norm, FFTW_RODFT00, FFTW_REDFT00); 01992 } 01993 01994 /********************************************************************/ 01995 01996 template <class SrcTraverser, class SrcAccessor, 01997 class DestTraverser, class DestAccessor> 01998 inline void 01999 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02000 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 02001 { 02002 fourierTransformRealEO(src.first, src.second, src.third, 02003 dest.first, dest.second, norm); 02004 } 02005 02006 template <class SrcTraverser, class SrcAccessor, 02007 class DestTraverser, class DestAccessor> 02008 inline void 02009 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02010 DestTraverser dul, DestAccessor dest, fftw_real norm) 02011 { 02012 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02013 norm, FFTW_REDFT00, FFTW_RODFT00); 02014 } 02015 02016 template <class DestTraverser, class DestAccessor> 02017 inline void 02018 fourierTransformRealEO( 02019 FFTWRealImage::const_traverser sul, 02020 FFTWRealImage::const_traverser slr, 02021 FFTWRealImage::Accessor src, 02022 DestTraverser dul, DestAccessor dest, fftw_real norm) 02023 { 02024 int w = slr.x - sul.x; 02025 02026 // test for right memory layout (fftw expects a width*height fftw_real array) 02027 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 02028 fourierTransformRealImpl(sul, slr, dul, dest, 02029 norm, FFTW_REDFT00, FFTW_RODFT00); 02030 else 02031 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02032 norm, FFTW_REDFT00, FFTW_RODFT00); 02033 } 02034 02035 /********************************************************************/ 02036 02037 template <class SrcTraverser, class SrcAccessor, 02038 class DestTraverser, class DestAccessor> 02039 inline void 02040 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02041 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 02042 { 02043 fourierTransformRealOO(src.first, src.second, src.third, 02044 dest.first, dest.second, norm); 02045 } 02046 02047 template <class SrcTraverser, class SrcAccessor, 02048 class DestTraverser, class DestAccessor> 02049 inline void 02050 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02051 DestTraverser dul, DestAccessor dest, fftw_real norm) 02052 { 02053 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02054 norm, FFTW_RODFT00, FFTW_RODFT00); 02055 } 02056 02057 template <class DestTraverser, class DestAccessor> 02058 inline void 02059 fourierTransformRealOO( 02060 FFTWRealImage::const_traverser sul, 02061 FFTWRealImage::const_traverser slr, 02062 FFTWRealImage::Accessor src, 02063 DestTraverser dul, DestAccessor dest, fftw_real norm) 02064 { 02065 int w = slr.x - sul.x; 02066 02067 // test for right memory layout (fftw expects a width*height fftw_real array) 02068 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 02069 fourierTransformRealImpl(sul, slr, dul, dest, 02070 norm, FFTW_RODFT00, FFTW_RODFT00); 02071 else 02072 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02073 norm, FFTW_RODFT00, FFTW_RODFT00); 02074 } 02075 02076 /*******************************************************************/ 02077 02078 template <class SrcTraverser, class SrcAccessor, 02079 class DestTraverser, class DestAccessor> 02080 void 02081 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02082 DestTraverser dul, DestAccessor dest, 02083 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy) 02084 { 02085 FFTWRealImage workImage(slr - sul); 02086 copyImage(srcIterRange(sul, slr, src), destImage(workImage)); 02087 FFTWRealImage const & cworkImage = workImage; 02088 fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), 02089 dul, dest, norm, kindx, kindy); 02090 } 02091 02092 02093 template <class DestTraverser, class DestAccessor> 02094 void 02095 fourierTransformRealImpl( 02096 FFTWRealImage::const_traverser sul, 02097 FFTWRealImage::const_traverser slr, 02098 DestTraverser dul, DestAccessor dest, 02099 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy) 02100 { 02101 int w = slr.x - sul.x; 02102 int h = slr.y - sul.y; 02103 BasicImage<fftw_real> res(w, h); 02104 02105 fftw_plan plan = fftw_plan_r2r_2d(h, w, 02106 (fftw_real *)&(*sul), (fftw_real *)res.begin(), 02107 kindy, kindx, FFTW_ESTIMATE); 02108 fftw_execute(plan); 02109 fftw_destroy_plan(plan); 02110 02111 if(norm != 1.0) 02112 transformImage(srcImageRange(res), destIter(dul, dest), 02113 std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm)); 02114 else 02115 copyImage(srcImageRange(res), destIter(dul, dest)); 02116 } 02117 02118 02119 //@} 02120 02121 } // namespace vigra 02122 02123 #endif // VIGRA_FFTW3_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|