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