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