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