[ 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.4.0, Dec 21 2005 ) */ 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 "vigra/stdimage.hxx" 00044 #include "vigra/copyimage.hxx" 00045 #include "vigra/transformimage.hxx" 00046 #include "vigra/combineimages.hxx" 00047 #include "vigra/numerictraits.hxx" 00048 #include "vigra/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 /********************************************************/ 01137 /* */ 01138 /* applyFourierFilter */ 01139 /* */ 01140 /********************************************************/ 01141 01142 /** \brief Apply a filter (defined in the frequency domain) to an image. 01143 01144 After transferring the image into the frequency domain, it is 01145 multiplied pixel-wise with the filter and transformed back. The 01146 result is always put into the given FFTWComplexImage 01147 <TT>destImg</TT> which must have the right size. Finally, the 01148 result will be normalized to compensate for the two FFTs. 01149 01150 The input and filter images can be scalar images (more precisely, 01151 <TT>SrcAccessor::value_type</TT> must be scalar) or 01152 FFTWComplexImages (in this case, one conversion is saved). 01153 01154 The DC entry of the filter must be in the upper left, which is the 01155 position where FFTW expects it (see \ref moveDCToUpperLeft()). 01156 01157 You can optionally pass two fftwnd_plans as last parameters, the 01158 forward plan and the in-place backward plan. Both must have been 01159 created for the right image size (see sample code). 01160 01161 <b> Declarations:</b> 01162 01163 pass arguments explicitly: 01164 \code 01165 namespace vigra { 01166 template <class SrcImageIterator, class SrcAccessor, 01167 class FilterImageIterator, class FilterAccessor> 01168 void applyFourierFilter(SrcImageIterator srcUpperLeft, 01169 SrcImageIterator srcLowerRight, SrcAccessor sa, 01170 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01171 FFTWComplexImage & destImg) 01172 01173 template <class SrcImageIterator, class SrcAccessor, 01174 class FilterImageIterator, class FilterAccessor> 01175 void applyFourierFilter(SrcImageIterator srcUpperLeft, 01176 SrcImageIterator srcLowerRight, SrcAccessor sa, 01177 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01178 FFTWComplexImage & destImg, 01179 const fftwnd_plan & forwardPlan, const fftwnd_plan & backwardPlan) 01180 } 01181 \endcode 01182 01183 use argument objects in conjunction with \ref ArgumentObjectFactories: 01184 \code 01185 namespace vigra { 01186 template <class SrcImageIterator, class SrcAccessor, 01187 class FilterImageIterator, class FilterAccessor> 01188 inline 01189 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01190 pair<FilterImageIterator, FilterAccessor> filter, 01191 FFTWComplexImage &destImg) 01192 01193 template <class SrcImageIterator, class SrcAccessor, 01194 class FilterImageIterator, class FilterAccessor> 01195 inline 01196 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01197 pair<FilterImageIterator, FilterAccessor> filter, 01198 FFTWComplexImage &destImg, 01199 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 01200 } 01201 \endcode 01202 01203 <b> Usage:</b> 01204 01205 <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br> 01206 Namespace: vigra 01207 01208 \code 01209 // create a Gaussian filter in Fourier space 01210 vigra::FImage gaussFilter(w, h), filter(w, h); 01211 for(int y=0; y<h; ++y) 01212 for(int x=0; x<w; ++x) 01213 { 01214 xx = float(x - w / 2) / w; 01215 yy = float(y - h / 2) / h; 01216 01217 gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale); 01218 } 01219 01220 // applyFourierFilter() expects the filter's DC in the upper left 01221 moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter)); 01222 01223 vigra::FFTWComplexImage result(w, h); 01224 01225 vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result); 01226 \endcode 01227 01228 For inspection of the result, \ref FFTWMagnitudeAccessor might be 01229 useful. If you want to apply the same filter repeatedly, it may be more 01230 efficient to use the FFTW functions directly with FFTW plans optimized 01231 for good performance. 01232 */ 01233 template <class SrcImageIterator, class SrcAccessor, 01234 class FilterImageIterator, class FilterAccessor, 01235 class DestImageIterator, class DestAccessor> 01236 inline 01237 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01238 pair<FilterImageIterator, FilterAccessor> filter, 01239 pair<DestImageIterator, DestAccessor> dest) 01240 { 01241 applyFourierFilter(src.first, src.second, src.third, 01242 filter.first, filter.second, 01243 dest.first, dest.second); 01244 } 01245 01246 template <class SrcImageIterator, class SrcAccessor, 01247 class FilterImageIterator, class FilterAccessor, 01248 class DestImageIterator, class DestAccessor> 01249 void applyFourierFilter(SrcImageIterator srcUpperLeft, 01250 SrcImageIterator srcLowerRight, SrcAccessor sa, 01251 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01252 DestImageIterator destUpperLeft, DestAccessor da) 01253 { 01254 // copy real input images into a complex one... 01255 int w= srcLowerRight.x - srcUpperLeft.x; 01256 int h= srcLowerRight.y - srcUpperLeft.y; 01257 01258 FFTWComplexImage workImage(w, h); 01259 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01260 destImage(workImage, FFTWWriteRealAccessor())); 01261 01262 // ...and call the impl 01263 FFTWComplexImage const & cworkImage = workImage; 01264 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01265 filterUpperLeft, fa, 01266 destUpperLeft, da); 01267 } 01268 01269 template <class FilterImageIterator, class FilterAccessor, 01270 class DestImageIterator, class DestAccessor> 01271 inline 01272 void applyFourierFilter( 01273 FFTWComplexImage::const_traverser srcUpperLeft, 01274 FFTWComplexImage::const_traverser srcLowerRight, 01275 FFTWComplexImage::ConstAccessor sa, 01276 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01277 DestImageIterator destUpperLeft, DestAccessor da) 01278 { 01279 int w = srcLowerRight.x - srcUpperLeft.x; 01280 int h = srcLowerRight.y - srcUpperLeft.y; 01281 01282 // test for right memory layout (fftw expects a 2*width*height floats array) 01283 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 01284 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa, 01285 filterUpperLeft, fa, 01286 destUpperLeft, da); 01287 else 01288 { 01289 FFTWComplexImage workImage(w, h); 01290 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01291 destImage(workImage)); 01292 01293 FFTWComplexImage const & cworkImage = workImage; 01294 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01295 filterUpperLeft, fa, 01296 destUpperLeft, da); 01297 } 01298 } 01299 01300 template <class FilterImageIterator, class FilterAccessor, 01301 class DestImageIterator, class DestAccessor> 01302 void applyFourierFilterImpl( 01303 FFTWComplexImage::const_traverser srcUpperLeft, 01304 FFTWComplexImage::const_traverser srcLowerRight, 01305 FFTWComplexImage::ConstAccessor sa, 01306 FilterImageIterator filterUpperLeft, FilterAccessor fa, 01307 DestImageIterator destUpperLeft, DestAccessor da) 01308 { 01309 int w = srcLowerRight.x - srcUpperLeft.x; 01310 int h = srcLowerRight.y - srcUpperLeft.y; 01311 01312 FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft); 01313 01314 // FFT from srcImage to complexResultImg 01315 fftw_plan forwardPlan= 01316 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft), 01317 (fftw_complex *)complexResultImg.begin(), 01318 FFTW_FORWARD, FFTW_ESTIMATE ); 01319 fftw_execute(forwardPlan); 01320 fftw_destroy_plan(forwardPlan); 01321 01322 // convolve in freq. domain (in complexResultImg) 01323 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa), 01324 destImage(complexResultImg), std::multiplies<FFTWComplex>()); 01325 01326 // FFT back into spatial domain (inplace in complexResultImg) 01327 fftw_plan backwardPlan= 01328 fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(), 01329 (fftw_complex *)complexResultImg.begin(), 01330 FFTW_BACKWARD, FFTW_ESTIMATE); 01331 fftw_execute(backwardPlan); 01332 fftw_destroy_plan(backwardPlan); 01333 01334 typedef typename 01335 NumericTraits<typename DestAccessor::value_type>::isScalar 01336 isScalarResult; 01337 01338 // normalization (after FFTs), maybe stripping imaginary part 01339 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da, 01340 isScalarResult()); 01341 } 01342 01343 template <class DestImageIterator, class DestAccessor> 01344 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage, 01345 DestImageIterator destUpperLeft, 01346 DestAccessor da, 01347 VigraFalseType) 01348 { 01349 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 01350 01351 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 01352 { 01353 DestImageIterator dIt= destUpperLeft; 01354 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 01355 { 01356 da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0); 01357 da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1); 01358 } 01359 } 01360 } 01361 01362 inline 01363 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 01364 FFTWComplexImage::traverser destUpperLeft, 01365 FFTWComplexImage::Accessor da, 01366 VigraFalseType) 01367 { 01368 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da), 01369 linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height()))); 01370 } 01371 01372 template <class DestImageIterator, class DestAccessor> 01373 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 01374 DestImageIterator destUpperLeft, 01375 DestAccessor da, 01376 VigraTrueType) 01377 { 01378 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 01379 01380 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 01381 { 01382 DestImageIterator dIt= destUpperLeft; 01383 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 01384 da.set(srcImage(x, y).re()*normFactor, dIt); 01385 } 01386 } 01387 01388 /**********************************************************/ 01389 /* */ 01390 /* applyFourierFilterFamily */ 01391 /* */ 01392 /**********************************************************/ 01393 01394 /** \brief Apply an array of filters (defined in the frequency domain) to an image. 01395 01396 This provides the same functionality as \ref applyFourierFilter(), 01397 but applying several filters at once allows to avoid 01398 repeated Fourier transforms of the source image. 01399 01400 Filters and result images must be stored in \ref vigra::ImageArray data 01401 structures. In contrast to \ref applyFourierFilter(), this function adjusts 01402 the size of the result images and the the length of the array. 01403 01404 <b> Declarations:</b> 01405 01406 pass arguments explicitly: 01407 \code 01408 namespace vigra { 01409 template <class SrcImageIterator, class SrcAccessor, class FilterType> 01410 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 01411 SrcImageIterator srcLowerRight, SrcAccessor sa, 01412 const ImageArray<FilterType> &filters, 01413 ImageArray<FFTWComplexImage> &results) 01414 } 01415 \endcode 01416 01417 use argument objects in conjunction with \ref ArgumentObjectFactories: 01418 \code 01419 namespace vigra { 01420 template <class SrcImageIterator, class SrcAccessor, class FilterType> 01421 inline 01422 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01423 const ImageArray<FilterType> &filters, 01424 ImageArray<FFTWComplexImage> &results) 01425 } 01426 \endcode 01427 01428 <b> Usage:</b> 01429 01430 <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br> 01431 Namespace: vigra 01432 01433 \code 01434 // assuming the presence of a real-valued image named "image" to 01435 // be filtered in this example 01436 01437 vigra::ImageArray<vigra::FImage> filters(16, image.size()); 01438 01439 for (int i=0; i<filters.size(); i++) 01440 // create some meaningful filters here 01441 createMyFilterOfScale(i, destImage(filters[i])); 01442 01443 vigra::ImageArray<vigra::FFTWComplexImage> results(); 01444 01445 vigra::applyFourierFilterFamily(srcImageRange(image), filters, results); 01446 \endcode 01447 */ 01448 template <class SrcImageIterator, class SrcAccessor, 01449 class FilterType, class DestImage> 01450 inline 01451 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01452 const ImageArray<FilterType> &filters, 01453 ImageArray<DestImage> &results) 01454 { 01455 applyFourierFilterFamily(src.first, src.second, src.third, 01456 filters, results); 01457 } 01458 01459 template <class SrcImageIterator, class SrcAccessor, 01460 class FilterType, class DestImage> 01461 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 01462 SrcImageIterator srcLowerRight, SrcAccessor sa, 01463 const ImageArray<FilterType> &filters, 01464 ImageArray<DestImage> &results) 01465 { 01466 int w= srcLowerRight.x - srcUpperLeft.x; 01467 int h= srcLowerRight.y - srcUpperLeft.y; 01468 01469 FFTWComplexImage workImage(w, h); 01470 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01471 destImage(workImage, FFTWWriteRealAccessor())); 01472 01473 FFTWComplexImage const & cworkImage = workImage; 01474 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01475 filters, results); 01476 } 01477 01478 template <class FilterType, class DestImage> 01479 inline 01480 void applyFourierFilterFamily( 01481 FFTWComplexImage::const_traverser srcUpperLeft, 01482 FFTWComplexImage::const_traverser srcLowerRight, 01483 FFTWComplexImage::ConstAccessor sa, 01484 const ImageArray<FilterType> &filters, 01485 ImageArray<DestImage> &results) 01486 { 01487 int w= srcLowerRight.x - srcUpperLeft.x; 01488 01489 // test for right memory layout (fftw expects a 2*width*height floats array) 01490 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 01491 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa, 01492 filters, results); 01493 else 01494 { 01495 int h = srcLowerRight.y - srcUpperLeft.y; 01496 FFTWComplexImage workImage(w, h); 01497 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01498 destImage(workImage)); 01499 01500 FFTWComplexImage const & cworkImage = workImage; 01501 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01502 filters, results); 01503 } 01504 } 01505 01506 template <class FilterType, class DestImage> 01507 void applyFourierFilterFamilyImpl( 01508 FFTWComplexImage::const_traverser srcUpperLeft, 01509 FFTWComplexImage::const_traverser srcLowerRight, 01510 FFTWComplexImage::ConstAccessor sa, 01511 const ImageArray<FilterType> &filters, 01512 ImageArray<DestImage> &results) 01513 { 01514 // make sure the filter images have the right dimensions 01515 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(), 01516 "applyFourierFilterFamily called with src image size != filters.imageSize()!"); 01517 01518 // make sure the result image array has the right dimensions 01519 results.resize(filters.size()); 01520 results.resizeImages(filters.imageSize()); 01521 01522 // FFT from srcImage to freqImage 01523 int w= srcLowerRight.x - srcUpperLeft.x; 01524 int h= srcLowerRight.y - srcUpperLeft.y; 01525 01526 FFTWComplexImage freqImage(w, h); 01527 FFTWComplexImage result(w, h); 01528 01529 fftw_plan forwardPlan= 01530 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft), 01531 (fftw_complex *)freqImage.begin(), 01532 FFTW_FORWARD, FFTW_ESTIMATE ); 01533 fftw_execute(forwardPlan); 01534 fftw_destroy_plan(forwardPlan); 01535 01536 fftw_plan backwardPlan= 01537 fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(), 01538 (fftw_complex *)result.begin(), 01539 FFTW_BACKWARD, FFTW_ESTIMATE ); 01540 typedef typename 01541 NumericTraits<typename DestImage::Accessor::value_type>::isScalar 01542 isScalarResult; 01543 01544 // convolve with filters in freq. domain 01545 for (unsigned int i= 0; i < filters.size(); i++) 01546 { 01547 combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]), 01548 destImage(result), std::multiplies<FFTWComplex>()); 01549 01550 // FFT back into spatial domain (inplace in destImage) 01551 fftw_execute(backwardPlan); 01552 01553 // normalization (after FFTs), maybe stripping imaginary part 01554 applyFourierFilterImplNormalization(result, 01555 results[i].upperLeft(), results[i].accessor(), 01556 isScalarResult()); 01557 } 01558 fftw_destroy_plan(backwardPlan); 01559 } 01560 01561 /********************************************************/ 01562 /* */ 01563 /* fourierTransformReal */ 01564 /* */ 01565 /********************************************************/ 01566 01567 /** \brief Real Fourier transforms for even and odd boundary conditions 01568 (aka. cosine and sine transforms). 01569 01570 01571 If the image is real and has even symmetry, its Fourier transform 01572 is also real and has even symmetry. The Fourier transform of a real image with odd 01573 symmetry is imaginary and has odd symmetry. In either case, only about a quarter 01574 of the pixels need to be stored because the rest can be calculated from the symmetry 01575 properties. This is especially useful, if the original image is implicitly assumed 01576 to have reflective or anti-reflective boundary conditions. Then the "negative" 01577 pixel locations are defined as 01578 01579 \code 01580 even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1) 01581 odd (anti-reflective boundary conditions): f[-1] = 0 01582 f[-x] = -f[x-2] (x = 2,...,N-1) 01583 \endcode 01584 01585 end similar at the other boundary (see the FFTW documentation for details). 01586 This has the advantage that more efficient Fourier transforms that use only 01587 real numbers can be implemented. These are also known as cosine and sine transforms 01588 respectively. 01589 01590 If you use the odd transform it is important to note that in the Fourier domain, 01591 the DC component is always zero and is therefore dropped from the data structure. 01592 This means that index 0 in an odd symmetric Fourier domain image refers to 01593 the <i>first</i> harmonic. This is especially important if an image is first 01594 cosine transformed (even symmetry), then in the Fourier domain multiplied 01595 with an odd symmetric filter (e.g. a first derivative) and finally transformed 01596 back to the spatial domain with a sine transform (odd symmetric). For this to work 01597 properly the image must be shifted left or up by one pixel (depending on whether 01598 the x- or y-axis is odd symmetric) before the inverse transform can be applied. 01599 (see example below). 01600 01601 The real Fourier transform functions are named <tt>fourierTransformReal??</tt> 01602 where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating 01603 whether the x- and y-axis is to be transformed using even or odd symmetry. 01604 The same functions can be used for both the forward and inverse transforms, 01605 only the normalization changes. For signal processing, the following 01606 normalization factors are most appropriate: 01607 01608 \code 01609 forward inverse 01610 ------------------------------------------------------------ 01611 X even, Y even 1.0 4.0 * (w-1) * (h-1) 01612 X even, Y odd -1.0 -4.0 * (w-1) * (h+1) 01613 X odd, Y even -1.0 -4.0 * (w+1) * (h-1) 01614 X odd, Y odd 1.0 4.0 * (w+1) * (h+1) 01615 \endcode 01616 01617 where <tt>w</tt> and <tt>h</tt> denote the image width and height. 01618 01619 <b> Declarations:</b> 01620 01621 pass arguments explicitly: 01622 \code 01623 namespace vigra { 01624 template <class SrcTraverser, class SrcAccessor, 01625 class DestTraverser, class DestAccessor> 01626 void 01627 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01628 DestTraverser dul, DestAccessor dest, fftw_real norm); 01629 01630 fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise 01631 } 01632 \endcode 01633 01634 01635 use argument objects in conjunction with \ref ArgumentObjectFactories: 01636 \code 01637 namespace vigra { 01638 template <class SrcTraverser, class SrcAccessor, 01639 class DestTraverser, class DestAccessor> 01640 void 01641 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01642 pair<DestTraverser, DestAccessor> dest, fftw_real norm); 01643 01644 fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise 01645 } 01646 \endcode 01647 01648 <b> Usage:</b> 01649 01650 <b>\#include</b> "<a href="fftw3_8hxx-source.html">vigra/fftw3.hxx</a>"<br> 01651 Namespace: vigra 01652 01653 \code 01654 vigra::FImage spatial(width,height), fourier(width,height); 01655 ... // fill image with data 01656 01657 // forward cosine transform == reflective boundary conditions 01658 fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0); 01659 01660 // multiply with a first derivative of Gaussian in x-direction 01661 for(int y = 0; y < height; ++y) 01662 { 01663 for(int x = 1; x < width; ++x) 01664 { 01665 double dx = x * M_PI / (width - 1); 01666 double dy = y * M_PI / (height - 1); 01667 fourier(x-1, y) = fourier(x, y) * dx * exp(-(dx*dx + dy*dy) * scale*scale / 2.0); 01668 } 01669 fourier(width-1, y) = 0.0; 01670 } 01671 01672 // inverse transform -- odd symmetry in x-direction, even in y, 01673 // due to symmetry of the filter 01674 fourierTransformRealOE(srcImageRange(fourier), destImage(spatial), 01675 (fftw_real)-4.0 * (width+1) * (height-1)); 01676 \endcode 01677 */ 01678 template <class SrcTraverser, class SrcAccessor, 01679 class DestTraverser, class DestAccessor> 01680 inline void 01681 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01682 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 01683 { 01684 fourierTransformRealEE(src.first, src.second, src.third, 01685 dest.first, dest.second, norm); 01686 } 01687 01688 template <class SrcTraverser, class SrcAccessor, 01689 class DestTraverser, class DestAccessor> 01690 inline void 01691 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01692 DestTraverser dul, DestAccessor dest, fftw_real norm) 01693 { 01694 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01695 norm, FFTW_REDFT00, FFTW_REDFT00); 01696 } 01697 01698 template <class DestTraverser, class DestAccessor> 01699 inline void 01700 fourierTransformRealEE( 01701 FFTWRealImage::const_traverser sul, 01702 FFTWRealImage::const_traverser slr, 01703 FFTWRealImage::Accessor src, 01704 DestTraverser dul, DestAccessor dest, fftw_real norm) 01705 { 01706 int w = slr.x - sul.x; 01707 01708 // test for right memory layout (fftw expects a width*height fftw_real array) 01709 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 01710 fourierTransformRealImpl(sul, slr, dul, dest, 01711 norm, FFTW_REDFT00, FFTW_REDFT00); 01712 else 01713 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01714 norm, FFTW_REDFT00, FFTW_REDFT00); 01715 } 01716 01717 /********************************************************************/ 01718 01719 template <class SrcTraverser, class SrcAccessor, 01720 class DestTraverser, class DestAccessor> 01721 inline void 01722 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01723 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 01724 { 01725 fourierTransformRealOE(src.first, src.second, src.third, 01726 dest.first, dest.second, norm); 01727 } 01728 01729 template <class SrcTraverser, class SrcAccessor, 01730 class DestTraverser, class DestAccessor> 01731 inline void 01732 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01733 DestTraverser dul, DestAccessor dest, fftw_real norm) 01734 { 01735 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01736 norm, FFTW_RODFT00, FFTW_REDFT00); 01737 } 01738 01739 template <class DestTraverser, class DestAccessor> 01740 inline void 01741 fourierTransformRealOE( 01742 FFTWRealImage::const_traverser sul, 01743 FFTWRealImage::const_traverser slr, 01744 FFTWRealImage::Accessor src, 01745 DestTraverser dul, DestAccessor dest, fftw_real norm) 01746 { 01747 int w = slr.x - sul.x; 01748 01749 // test for right memory layout (fftw expects a width*height fftw_real array) 01750 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 01751 fourierTransformRealImpl(sul, slr, dul, dest, 01752 norm, FFTW_RODFT00, FFTW_REDFT00); 01753 else 01754 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01755 norm, FFTW_RODFT00, FFTW_REDFT00); 01756 } 01757 01758 /********************************************************************/ 01759 01760 template <class SrcTraverser, class SrcAccessor, 01761 class DestTraverser, class DestAccessor> 01762 inline void 01763 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01764 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 01765 { 01766 fourierTransformRealEO(src.first, src.second, src.third, 01767 dest.first, dest.second, norm); 01768 } 01769 01770 template <class SrcTraverser, class SrcAccessor, 01771 class DestTraverser, class DestAccessor> 01772 inline void 01773 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01774 DestTraverser dul, DestAccessor dest, fftw_real norm) 01775 { 01776 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01777 norm, FFTW_REDFT00, FFTW_RODFT00); 01778 } 01779 01780 template <class DestTraverser, class DestAccessor> 01781 inline void 01782 fourierTransformRealEO( 01783 FFTWRealImage::const_traverser sul, 01784 FFTWRealImage::const_traverser slr, 01785 FFTWRealImage::Accessor src, 01786 DestTraverser dul, DestAccessor dest, fftw_real norm) 01787 { 01788 int w = slr.x - sul.x; 01789 01790 // test for right memory layout (fftw expects a width*height fftw_real array) 01791 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 01792 fourierTransformRealImpl(sul, slr, dul, dest, 01793 norm, FFTW_REDFT00, FFTW_RODFT00); 01794 else 01795 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01796 norm, FFTW_REDFT00, FFTW_RODFT00); 01797 } 01798 01799 /********************************************************************/ 01800 01801 template <class SrcTraverser, class SrcAccessor, 01802 class DestTraverser, class DestAccessor> 01803 inline void 01804 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 01805 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 01806 { 01807 fourierTransformRealOO(src.first, src.second, src.third, 01808 dest.first, dest.second, norm); 01809 } 01810 01811 template <class SrcTraverser, class SrcAccessor, 01812 class DestTraverser, class DestAccessor> 01813 inline void 01814 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01815 DestTraverser dul, DestAccessor dest, fftw_real norm) 01816 { 01817 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01818 norm, FFTW_RODFT00, FFTW_RODFT00); 01819 } 01820 01821 template <class DestTraverser, class DestAccessor> 01822 inline void 01823 fourierTransformRealOO( 01824 FFTWRealImage::const_traverser sul, 01825 FFTWRealImage::const_traverser slr, 01826 FFTWRealImage::Accessor src, 01827 DestTraverser dul, DestAccessor dest, fftw_real norm) 01828 { 01829 int w = slr.x - sul.x; 01830 01831 // test for right memory layout (fftw expects a width*height fftw_real array) 01832 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 01833 fourierTransformRealImpl(sul, slr, dul, dest, 01834 norm, FFTW_RODFT00, FFTW_RODFT00); 01835 else 01836 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 01837 norm, FFTW_RODFT00, FFTW_RODFT00); 01838 } 01839 01840 /*******************************************************************/ 01841 01842 template <class SrcTraverser, class SrcAccessor, 01843 class DestTraverser, class DestAccessor> 01844 void 01845 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 01846 DestTraverser dul, DestAccessor dest, 01847 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy) 01848 { 01849 FFTWRealImage workImage(slr - sul); 01850 copyImage(srcIterRange(sul, slr, src), destImage(workImage)); 01851 FFTWRealImage const & cworkImage = workImage; 01852 fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), 01853 dul, dest, norm, kindx, kindy); 01854 } 01855 01856 01857 template <class DestTraverser, class DestAccessor> 01858 void 01859 fourierTransformRealImpl( 01860 FFTWRealImage::const_traverser sul, 01861 FFTWRealImage::const_traverser slr, 01862 DestTraverser dul, DestAccessor dest, 01863 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy) 01864 { 01865 int w = slr.x - sul.x; 01866 int h = slr.y - sul.y; 01867 BasicImage<fftw_real> res(w, h); 01868 01869 fftw_plan plan = fftw_plan_r2r_2d(h, w, 01870 (fftw_real *)&(*sul), (fftw_real *)res.begin(), 01871 kindy, kindx, FFTW_ESTIMATE); 01872 fftw_execute(plan); 01873 fftw_destroy_plan(plan); 01874 01875 if(norm != 1.0) 01876 transformImage(srcImageRange(res), destIter(dul, dest), 01877 std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm)); 01878 else 01879 copyImage(srcImageRange(res), destIter(dul, dest)); 01880 } 01881 01882 01883 //@} 01884 01885 } // namespace vigra 01886 01887 #endif // VIGRA_FFTW3_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|