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