[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/fftw3.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2004 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_FFTW3_HXX 00037 #define VIGRA_FFTW3_HXX 00038 00039 #include <cmath> 00040 #include <functional> 00041 #include <complex> 00042 #include "stdimage.hxx" 00043 #include "copyimage.hxx" 00044 #include "transformimage.hxx" 00045 #include "combineimages.hxx" 00046 #include "numerictraits.hxx" 00047 #include "imagecontainer.hxx" 00048 #include <fftw3.h> 00049 00050 namespace vigra { 00051 00052 typedef double fftw_real; 00053 00054 template <class T> 00055 struct FFTWReal; 00056 00057 template <> 00058 struct FFTWReal<fftw_complex> 00059 { 00060 typedef double type; 00061 }; 00062 00063 template <> 00064 struct FFTWReal<fftwf_complex> 00065 { 00066 typedef float type; 00067 }; 00068 00069 template <> 00070 struct FFTWReal<fftwl_complex> 00071 { 00072 typedef long double type; 00073 }; 00074 00075 template <class T> 00076 struct FFTWReal2Complex; 00077 00078 template <> 00079 struct FFTWReal2Complex<double> 00080 { 00081 typedef fftw_complex type; 00082 typedef fftw_plan plan_type; 00083 }; 00084 00085 template <> 00086 struct FFTWReal2Complex<float> 00087 { 00088 typedef fftwf_complex type; 00089 typedef fftwf_plan plan_type; 00090 }; 00091 00092 template <> 00093 struct FFTWReal2Complex<long double> 00094 { 00095 typedef fftwl_complex type; 00096 typedef fftwl_plan plan_type; 00097 }; 00098 00099 /********************************************************/ 00100 /* */ 00101 /* FFTWComplex */ 00102 /* */ 00103 /********************************************************/ 00104 00105 /** \brief Wrapper class for the FFTW complex types '<TT>fftw_complex</TT>'. 00106 00107 This class encapsulates the low-level complex number types provided by the 00108 <a href="http://www.fftw.org/">FFTW Fast Fourier Transform</a> library (i.e. 00109 '<TT>fftw_complex</TT>', '<TT>fftwf_complex</TT>', '<TT>fftwl_complex</TT>'). 00110 In particular, it provides constructors, member functions and 00111 \ref FFTWComplexOperators "arithmetic operators" that make FFTW complex numbers 00112 compatible with <tt>std::complex</tt>. In addition, the class defines 00113 transformations to polar coordinates and \ref FFTWComplexAccessors "accessors". 00114 00115 FFTWComplex implements the concepts \ref AlgebraicField and 00116 \ref DivisionAlgebra. The standard image types <tt>FFTWRealImage</tt> 00117 and <tt>FFTWComplexImage</tt> are defined. 00118 00119 <b>See also:</b> 00120 <ul> 00121 <li> \ref FFTWComplexTraits 00122 <li> \ref FFTWComplexOperators 00123 <li> \ref FFTWComplexAccessors 00124 </ul> 00125 00126 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 00127 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 00128 Namespace: vigra 00129 */ 00130 template <class Real = double> 00131 class FFTWComplex 00132 { 00133 public: 00134 /** The wrapped complex type 00135 */ 00136 typedef typename FFTWReal2Complex<Real>::type complex_type; 00137 00138 /** The complex' component type, as defined in '<TT>fftw3.h</TT>' 00139 */ 00140 typedef Real value_type; 00141 00142 /** reference type (result of operator[]) 00143 */ 00144 typedef value_type & reference; 00145 00146 /** const reference type (result of operator[] const) 00147 */ 00148 typedef value_type const & const_reference; 00149 00150 /** iterator type (result of begin() ) 00151 */ 00152 typedef value_type * iterator; 00153 00154 /** const iterator type (result of begin() const) 00155 */ 00156 typedef value_type const * const_iterator; 00157 00158 /** The norm type (result of magnitude()) 00159 */ 00160 typedef value_type NormType; 00161 00162 /** The squared norm type (result of squaredMagnitde()) 00163 */ 00164 typedef value_type SquaredNormType; 00165 00166 /** Construct from real and imaginary part. 00167 Default: 0. 00168 */ 00169 FFTWComplex(value_type const & re = 0.0, value_type const & im = 0.0) 00170 { 00171 data_[0] = re; 00172 data_[1] = im; 00173 } 00174 00175 /** Copy constructor. 00176 */ 00177 FFTWComplex(FFTWComplex const & o) 00178 { 00179 data_[0] = o.data_[0]; 00180 data_[1] = o.data_[1]; 00181 } 00182 00183 /** Copy constructor. 00184 */ 00185 template <class U> 00186 FFTWComplex(FFTWComplex<U> const & o) 00187 { 00188 data_[0] = (Real)o.real(); 00189 data_[1] = (Real)o.imag(); 00190 } 00191 00192 /** Construct from plain <TT>fftw_complex</TT>. 00193 */ 00194 FFTWComplex(fftw_complex const & o) 00195 { 00196 data_[0] = (Real)o[0]; 00197 data_[1] = (Real)o[1]; 00198 } 00199 00200 /** Construct from plain <TT>fftwf_complex</TT>. 00201 */ 00202 FFTWComplex(fftwf_complex const & o) 00203 { 00204 data_[0] = (Real)o[0]; 00205 data_[1] = (Real)o[1]; 00206 } 00207 00208 /** Construct from plain <TT>fftwl_complex</TT>. 00209 */ 00210 FFTWComplex(fftwl_complex const & o) 00211 { 00212 data_[0] = (Real)o[0]; 00213 data_[1] = (Real)o[1]; 00214 } 00215 00216 /** Construct from std::complex. 00217 */ 00218 template <class T> 00219 FFTWComplex(std::complex<T> const & o) 00220 { 00221 data_[0] = (Real)o.real(); 00222 data_[1] = (Real)o.imag(); 00223 } 00224 00225 /** Construct from TinyVector. 00226 */ 00227 template <class T> 00228 FFTWComplex(TinyVector<T, 2> const & o) 00229 { 00230 data_[0] = (Real)o[0]; 00231 data_[1] = (Real)o[1]; 00232 } 00233 00234 /** Assignment. 00235 */ 00236 FFTWComplex& operator=(FFTWComplex const & o) 00237 { 00238 data_[0] = o.data_[0]; 00239 data_[1] = o.data_[1]; 00240 return *this; 00241 } 00242 00243 /** Assignment. 00244 */ 00245 template <class U> 00246 FFTWComplex& operator=(FFTWComplex<U> const & o) 00247 { 00248 data_[0] = (Real)o.real(); 00249 data_[1] = (Real)o.imag(); 00250 return *this; 00251 } 00252 00253 /** Assignment. 00254 */ 00255 FFTWComplex& operator=(fftw_complex const & o) 00256 { 00257 data_[0] = (Real)o[0]; 00258 data_[1] = (Real)o[1]; 00259 return *this; 00260 } 00261 00262 /** Assignment. 00263 */ 00264 FFTWComplex& operator=(fftwf_complex const & o) 00265 { 00266 data_[0] = (Real)o[0]; 00267 data_[1] = (Real)o[1]; 00268 return *this; 00269 } 00270 00271 /** Assignment. 00272 */ 00273 FFTWComplex& operator=(fftwl_complex const & o) 00274 { 00275 data_[0] = (Real)o[0]; 00276 data_[1] = (Real)o[1]; 00277 return *this; 00278 } 00279 00280 /** Assignment. 00281 */ 00282 FFTWComplex& operator=(double o) 00283 { 00284 data_[0] = (Real)o; 00285 data_[1] = 0.0; 00286 return *this; 00287 } 00288 00289 /** Assignment. 00290 */ 00291 FFTWComplex& operator=(float o) 00292 { 00293 data_[0] = (Real)o; 00294 data_[1] = 0.0; 00295 return *this; 00296 } 00297 00298 /** Assignment. 00299 */ 00300 FFTWComplex& operator=(long double o) 00301 { 00302 data_[0] = (Real)o; 00303 data_[1] = 0.0; 00304 return *this; 00305 } 00306 00307 /** Assignment. 00308 */ 00309 template <class T> 00310 FFTWComplex& operator=(TinyVector<T, 2> const & o) 00311 { 00312 data_[0] = (Real)o[0]; 00313 data_[1] = (Real)o[1]; 00314 return *this; 00315 } 00316 00317 /** Assignment. 00318 */ 00319 template <class T> 00320 FFTWComplex& operator=(std::complex<T> const & o) 00321 { 00322 data_[0] = (Real)o.real(); 00323 data_[1] = (Real)o.imag(); 00324 return *this; 00325 } 00326 00327 reference re() 00328 { return data_[0]; } 00329 00330 const_reference re() const 00331 { return data_[0]; } 00332 00333 reference real() 00334 { return data_[0]; } 00335 00336 const_reference real() const 00337 { return data_[0]; } 00338 00339 reference im() 00340 { return data_[1]; } 00341 00342 const_reference im() const 00343 { return data_[1]; } 00344 00345 reference imag() 00346 { return data_[1]; } 00347 00348 const_reference imag() const 00349 { return data_[1]; } 00350 00351 /** Unary negation. 00352 */ 00353 FFTWComplex operator-() const 00354 { return FFTWComplex(-data_[0], -data_[1]); } 00355 00356 /** Squared magnitude x*conj(x) 00357 */ 00358 SquaredNormType squaredMagnitude() const 00359 { return data_[0]*data_[0]+data_[1]*data_[1]; } 00360 00361 /** Magnitude (length of radius vector). 00362 */ 00363 NormType magnitude() const 00364 { return VIGRA_CSTD::sqrt(squaredMagnitude()); } 00365 00366 /** Phase angle. 00367 */ 00368 value_type phase() const 00369 { return VIGRA_CSTD::atan2(data_[1], data_[0]); } 00370 00371 /** Access components as if number were a vector. 00372 */ 00373 reference operator[](int i) 00374 { return data_[i]; } 00375 00376 /** Read components as if number were a vector. 00377 */ 00378 const_reference operator[](int i) const 00379 { return data_[i]; } 00380 00381 /** Length of complex number (always 2). 00382 */ 00383 int size() const 00384 { return 2; } 00385 00386 iterator begin() 00387 { return data_; } 00388 00389 iterator end() 00390 { return data_ + 2; } 00391 00392 const_iterator begin() const 00393 { return data_; } 00394 00395 const_iterator end() const 00396 { return data_ + 2; } 00397 00398 private: 00399 complex_type data_; 00400 }; 00401 00402 /********************************************************/ 00403 /* */ 00404 /* FFTWComplexTraits */ 00405 /* */ 00406 /********************************************************/ 00407 00408 /** \page FFTWComplexTraits Numeric and Promote Traits of FFTWComplex 00409 00410 The numeric and promote traits for fftw_complex and FFTWComplex follow 00411 the general specifications for \ref NumericPromotionTraits and 00412 \ref AlgebraicField. They are explicitly specialized for the types 00413 involved: 00414 00415 \code 00416 00417 template<> 00418 struct NumericTraits<fftw_complex> 00419 { 00420 typedef fftw_complex Promote; 00421 typedef fftw_complex RealPromote; 00422 typedef fftw_complex ComplexPromote; 00423 typedef fftw_real ValueType; 00424 00425 typedef VigraFalseType isIntegral; 00426 typedef VigraFalseType isScalar; 00427 typedef VigraFalseType isOrdered; 00428 typedef VigraTrueType isComplex; 00429 00430 // etc. 00431 }; 00432 00433 template<class Real> 00434 struct NumericTraits<FFTWComplex<Real> > 00435 { 00436 typedef FFTWComplex<Real> Promote; 00437 typedef FFTWComplex<Real> RealPromote; 00438 typedef FFTWComplex<Real> ComplexPromote; 00439 typedef Real ValueType; 00440 00441 typedef VigraFalseType isIntegral; 00442 typedef VigraFalseType isScalar; 00443 typedef VigraFalseType isOrdered; 00444 typedef VigraTrueType isComplex; 00445 00446 // etc. 00447 }; 00448 00449 template<> 00450 struct NormTraits<fftw_complex> 00451 { 00452 typedef fftw_complex Type; 00453 typedef fftw_real SquaredNormType; 00454 typedef fftw_real NormType; 00455 }; 00456 00457 template<class Real> 00458 struct NormTraits<FFTWComplex> 00459 { 00460 typedef FFTWComplex Type; 00461 typedef fftw_real SquaredNormType; 00462 typedef fftw_real NormType; 00463 }; 00464 00465 template <> 00466 struct PromoteTraits<fftw_complex, fftw_complex> 00467 { 00468 typedef fftw_complex Promote; 00469 }; 00470 00471 template <> 00472 struct PromoteTraits<fftw_complex, double> 00473 { 00474 typedef fftw_complex Promote; 00475 }; 00476 00477 template <> 00478 struct PromoteTraits<double, fftw_complex> 00479 { 00480 typedef fftw_complex Promote; 00481 }; 00482 00483 template <class Real> 00484 struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> > 00485 { 00486 typedef FFTWComplex<Real> Promote; 00487 }; 00488 00489 template <class Real> 00490 struct PromoteTraits<FFTWComplex<Real>, double> 00491 { 00492 typedef FFTWComplex<Real> Promote; 00493 }; 00494 00495 template <class Real> 00496 struct PromoteTraits<double, FFTWComplex<Real> > 00497 { 00498 typedef FFTWComplex<Real> Promote; 00499 }; 00500 \endcode 00501 00502 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 00503 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 00504 Namespace: vigra 00505 00506 */ 00507 template<> 00508 struct NumericTraits<fftw_complex> 00509 { 00510 typedef fftw_complex Type; 00511 typedef fftw_complex Promote; 00512 typedef fftw_complex RealPromote; 00513 typedef fftw_complex ComplexPromote; 00514 typedef fftw_real ValueType; 00515 00516 typedef VigraFalseType isIntegral; 00517 typedef VigraFalseType isScalar; 00518 typedef NumericTraits<fftw_real>::isSigned isSigned; 00519 typedef VigraFalseType isOrdered; 00520 typedef VigraTrueType isComplex; 00521 00522 static FFTWComplex<> zero() { return FFTWComplex<>(0.0, 0.0); } 00523 static FFTWComplex<> one() { return FFTWComplex<>(1.0, 0.0); } 00524 static FFTWComplex<> nonZero() { return one(); } 00525 00526 static const Promote & toPromote(const Type & v) { return v; } 00527 static const RealPromote & toRealPromote(const Type & v) { return v; } 00528 static const Type & fromPromote(const Promote & v) { return v; } 00529 static const Type & fromRealPromote(const RealPromote & v) { return v; } 00530 }; 00531 00532 template<class Real> 00533 struct NumericTraits<FFTWComplex<Real> > 00534 { 00535 typedef FFTWComplex<Real> Type; 00536 typedef FFTWComplex<Real> Promote; 00537 typedef FFTWComplex<Real> RealPromote; 00538 typedef FFTWComplex<Real> ComplexPromote; 00539 typedef typename Type::value_type ValueType; 00540 00541 typedef VigraFalseType isIntegral; 00542 typedef VigraFalseType isScalar; 00543 typedef typename NumericTraits<ValueType>::isSigned isSigned; 00544 typedef VigraFalseType isOrdered; 00545 typedef VigraTrueType isComplex; 00546 00547 static Type zero() { return Type(0.0, 0.0); } 00548 static Type one() { return Type(1.0, 0.0); } 00549 static Type nonZero() { return one(); } 00550 00551 static const Promote & toPromote(const Type & v) { return v; } 00552 static const RealPromote & toRealPromote(const Type & v) { return v; } 00553 static const Type & fromPromote(const Promote & v) { return v; } 00554 static const Type & fromRealPromote(const RealPromote & v) { return v; } 00555 }; 00556 00557 template<> 00558 struct NormTraits<fftw_complex> 00559 { 00560 typedef fftw_complex Type; 00561 typedef fftw_real SquaredNormType; 00562 typedef fftw_real NormType; 00563 }; 00564 00565 template<class Real> 00566 struct NormTraits<FFTWComplex<Real> > 00567 { 00568 typedef FFTWComplex<Real> Type; 00569 typedef typename Type::SquaredNormType SquaredNormType; 00570 typedef typename Type::NormType NormType; 00571 }; 00572 00573 template <> 00574 struct PromoteTraits<fftw_complex, fftw_complex> 00575 { 00576 typedef fftw_complex Promote; 00577 }; 00578 00579 template <> 00580 struct PromoteTraits<fftw_complex, double> 00581 { 00582 typedef fftw_complex Promote; 00583 }; 00584 00585 template <> 00586 struct PromoteTraits<double, fftw_complex> 00587 { 00588 typedef fftw_complex Promote; 00589 }; 00590 00591 template <class Real> 00592 struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> > 00593 { 00594 typedef FFTWComplex<Real> Promote; 00595 }; 00596 00597 template <class Real> 00598 struct PromoteTraits<FFTWComplex<Real>, double> 00599 { 00600 typedef FFTWComplex<Real> Promote; 00601 }; 00602 00603 template <class Real> 00604 struct PromoteTraits<double, FFTWComplex<Real> > 00605 { 00606 typedef FFTWComplex<Real> Promote; 00607 }; 00608 00609 template<class T> 00610 struct CanSkipInitialization<std::complex<T> > 00611 { 00612 typedef typename CanSkipInitialization<T>::type type; 00613 static const bool value = type::asBool; 00614 }; 00615 00616 template<class Real> 00617 struct CanSkipInitialization<FFTWComplex<Real> > 00618 { 00619 typedef typename CanSkipInitialization<Real>::type type; 00620 static const bool value = type::asBool; 00621 }; 00622 00623 namespace multi_math { 00624 00625 template <class ARG> 00626 struct MultiMathOperand; 00627 00628 template <class Real> 00629 struct MultiMathOperand<FFTWComplex<Real> > 00630 { 00631 typedef MultiMathOperand<FFTWComplex<Real> > AllowOverload; 00632 typedef FFTWComplex<Real> result_type; 00633 00634 static const int ndim = 0; 00635 00636 MultiMathOperand(FFTWComplex<Real> const & v) 00637 : v_(v) 00638 {} 00639 00640 template <class SHAPE> 00641 bool checkShape(SHAPE const &) const 00642 { 00643 return true; 00644 } 00645 00646 template <class SHAPE> 00647 FFTWComplex<Real> const & operator[](SHAPE const &) const 00648 { 00649 return v_; 00650 } 00651 00652 void inc(unsigned int /*LEVEL*/) const 00653 {} 00654 00655 void reset(unsigned int /*LEVEL*/) const 00656 {} 00657 00658 FFTWComplex<Real> const & operator*() const 00659 { 00660 return v_; 00661 } 00662 00663 FFTWComplex<Real> v_; 00664 }; 00665 00666 } // namespace multi_math 00667 00668 template<class Ty> 00669 class FFTWAllocator 00670 { 00671 public: 00672 typedef std::size_t size_type; 00673 typedef std::ptrdiff_t difference_type; 00674 typedef Ty *pointer; 00675 typedef const Ty *const_pointer; 00676 typedef Ty& reference; 00677 typedef const Ty& const_reference; 00678 typedef Ty value_type; 00679 00680 pointer address(reference val) const 00681 { return &val; } 00682 00683 const_pointer address(const_reference val) const 00684 { return &val; } 00685 00686 template<class Other> 00687 struct rebind 00688 { 00689 typedef FFTWAllocator<Other> other; 00690 }; 00691 00692 FFTWAllocator() throw() 00693 {} 00694 00695 template<class Other> 00696 FFTWAllocator(const FFTWAllocator<Other>& right) throw() 00697 {} 00698 00699 template<class Other> 00700 FFTWAllocator& operator=(const FFTWAllocator<Other>& right) 00701 { 00702 return *this; 00703 } 00704 00705 pointer allocate(size_type count, void * = 0) 00706 { 00707 return (pointer)fftw_malloc(count * sizeof(Ty)); 00708 } 00709 00710 void deallocate(pointer ptr, size_type count) 00711 { 00712 fftw_free(ptr); 00713 } 00714 00715 void construct(pointer ptr, const Ty& val) 00716 { 00717 new(ptr) Ty(val); 00718 00719 } 00720 00721 void destroy(pointer ptr) 00722 { 00723 ptr->~Ty(); 00724 } 00725 00726 size_type max_size() const throw() 00727 { 00728 return NumericTraits<std::ptrdiff_t>::max() / sizeof(Ty); 00729 } 00730 }; 00731 00732 } // namespace vigra 00733 00734 namespace std { 00735 00736 template<class Real> 00737 class allocator<vigra::FFTWComplex<Real> > 00738 { 00739 public: 00740 typedef vigra::FFTWComplex<Real> value_type; 00741 typedef size_t size_type; 00742 typedef std::ptrdiff_t difference_type; 00743 typedef value_type *pointer; 00744 typedef const value_type *const_pointer; 00745 typedef value_type& reference; 00746 typedef const value_type& const_reference; 00747 00748 pointer address(reference val) const 00749 { return &val; } 00750 00751 const_pointer address(const_reference val) const 00752 { return &val; } 00753 00754 template<class Other> 00755 struct rebind 00756 { 00757 typedef allocator<Other> other; 00758 }; 00759 00760 allocator() throw() 00761 {} 00762 00763 template<class Other> 00764 allocator(const allocator<Other>& right) throw() 00765 {} 00766 00767 template<class Other> 00768 allocator& operator=(const allocator<Other>& right) 00769 { 00770 return *this; 00771 } 00772 00773 pointer allocate(size_type count, void * = 0) 00774 { 00775 return (pointer)fftw_malloc(count * sizeof(value_type)); 00776 } 00777 00778 void deallocate(pointer ptr, size_type /*count*/) 00779 { 00780 fftw_free(ptr); 00781 } 00782 00783 void construct(pointer ptr, const value_type& val) 00784 { 00785 new(ptr) value_type(val); 00786 00787 } 00788 00789 void destroy(pointer ptr) 00790 { 00791 ptr->~value_type(); 00792 } 00793 00794 size_type max_size() const throw() 00795 { 00796 return vigra::NumericTraits<std::ptrdiff_t>::max() / sizeof(value_type); 00797 } 00798 }; 00799 00800 } // namespace std 00801 00802 namespace vigra { 00803 00804 /********************************************************/ 00805 /* */ 00806 /* FFTWComplex Operations */ 00807 /* */ 00808 /********************************************************/ 00809 00810 /** \addtogroup FFTWComplexOperators Functions for FFTWComplex 00811 00812 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 00813 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 00814 00815 These functions fulfill the requirements of an Algebraic Field. 00816 Return types are determined according to \ref FFTWComplexTraits. 00817 00818 Namespace: vigra 00819 <p> 00820 00821 */ 00822 //@{ 00823 /// equal 00824 template <class R> 00825 inline bool operator ==(FFTWComplex<R> const &a, const FFTWComplex<R> &b) { 00826 return a.re() == b.re() && a.im() == b.im(); 00827 } 00828 00829 template <class R> 00830 inline bool operator ==(FFTWComplex<R> const &a, double b) { 00831 return a.re() == b && a.im() == 0.0; 00832 } 00833 00834 template <class R> 00835 inline bool operator ==(double a, const FFTWComplex<R> &b) { 00836 return a == b.re() && 0.0 == b.im(); 00837 } 00838 00839 /// not equal 00840 template <class R> 00841 inline bool operator !=(FFTWComplex<R> const &a, const FFTWComplex<R> &b) { 00842 return a.re() != b.re() || a.im() != b.im(); 00843 } 00844 00845 /// not equal 00846 template <class R> 00847 inline bool operator !=(FFTWComplex<R> const &a, double b) { 00848 return a.re() != b || a.im() != 0.0; 00849 } 00850 00851 /// not equal 00852 template <class R> 00853 inline bool operator !=(double a, const FFTWComplex<R> &b) { 00854 return a != b.re() || 0.0 != b.im(); 00855 } 00856 00857 /// add-assignment 00858 template <class R> 00859 inline FFTWComplex<R> & operator +=(FFTWComplex<R> &a, const FFTWComplex<R> &b) { 00860 a.re() += b.re(); 00861 a.im() += b.im(); 00862 return a; 00863 } 00864 00865 /// subtract-assignment 00866 template <class R> 00867 inline FFTWComplex<R> & operator -=(FFTWComplex<R> &a, const FFTWComplex<R> &b) { 00868 a.re() -= b.re(); 00869 a.im() -= b.im(); 00870 return a; 00871 } 00872 00873 /// multiply-assignment 00874 template <class R> 00875 inline FFTWComplex<R> & operator *=(FFTWComplex<R> &a, const FFTWComplex<R> &b) { 00876 typename FFTWComplex<R>::value_type t = a.re()*b.re()-a.im()*b.im(); 00877 a.im() = a.re()*b.im()+a.im()*b.re(); 00878 a.re() = t; 00879 return a; 00880 } 00881 00882 /// divide-assignment 00883 template <class R> 00884 inline FFTWComplex<R> & operator /=(FFTWComplex<R> &a, const FFTWComplex<R> &b) { 00885 typename FFTWComplex<R>::value_type sm = b.squaredMagnitude(); 00886 typename FFTWComplex<R>::value_type t = (a.re()*b.re()+a.im()*b.im())/sm; 00887 a.im() = (b.re()*a.im()-a.re()*b.im())/sm; 00888 a.re() = t; 00889 return a; 00890 } 00891 00892 /// add-assignment with scalar double 00893 template <class R> 00894 inline FFTWComplex<R> & operator +=(FFTWComplex<R> &a, double b) { 00895 a.re() += (R)b; 00896 return a; 00897 } 00898 00899 /// subtract-assignment with scalar double 00900 template <class R> 00901 inline FFTWComplex<R> & operator -=(FFTWComplex<R> &a, double b) { 00902 a.re() -= (R)b; 00903 return a; 00904 } 00905 00906 /// multiply-assignment with scalar double 00907 template <class R> 00908 inline FFTWComplex<R> & operator *=(FFTWComplex<R> &a, double b) { 00909 a.re() *= (R)b; 00910 a.im() *= (R)b; 00911 return a; 00912 } 00913 00914 /// divide-assignment with scalar double 00915 template <class R> 00916 inline FFTWComplex<R> & operator /=(FFTWComplex<R> &a, double b) { 00917 a.re() /= (R)b; 00918 a.im() /= (R)b; 00919 return a; 00920 } 00921 00922 /// addition 00923 template <class R> 00924 inline FFTWComplex<R> operator +(FFTWComplex<R> a, const FFTWComplex<R> &b) { 00925 a += b; 00926 return a; 00927 } 00928 00929 /// right addition with scalar double 00930 template <class R> 00931 inline FFTWComplex<R> operator +(FFTWComplex<R> a, double b) { 00932 a += b; 00933 return a; 00934 } 00935 00936 /// left addition with scalar double 00937 template <class R> 00938 inline FFTWComplex<R> operator +(double a, FFTWComplex<R> b) { 00939 b += a; 00940 return b; 00941 } 00942 00943 /// subtraction 00944 template <class R> 00945 inline FFTWComplex<R> operator -(FFTWComplex<R> a, const FFTWComplex<R> &b) { 00946 a -= b; 00947 return a; 00948 } 00949 00950 /// right subtraction with scalar double 00951 template <class R> 00952 inline FFTWComplex<R> operator -(FFTWComplex<R> a, double b) { 00953 a -= b; 00954 return a; 00955 } 00956 00957 /// left subtraction with scalar double 00958 template <class R> 00959 inline FFTWComplex<R> operator -(double a, FFTWComplex<R> const & b) { 00960 return (-b) += a; 00961 } 00962 00963 /// multiplication 00964 template <class R> 00965 inline FFTWComplex<R> operator *(FFTWComplex<R> a, const FFTWComplex<R> &b) { 00966 a *= b; 00967 return a; 00968 } 00969 00970 /// right multiplication with scalar double 00971 template <class R> 00972 inline FFTWComplex<R> operator *(FFTWComplex<R> a, double b) { 00973 a *= b; 00974 return a; 00975 } 00976 00977 /// left multiplication with scalar double 00978 template <class R> 00979 inline FFTWComplex<R> operator *(double a, FFTWComplex<R> b) { 00980 b *= a; 00981 return b; 00982 } 00983 00984 /// division 00985 template <class R> 00986 inline FFTWComplex<R> operator /(FFTWComplex<R> a, const FFTWComplex<R> &b) { 00987 a /= b; 00988 return a; 00989 } 00990 00991 /// right division with scalar double 00992 template <class R> 00993 inline FFTWComplex<R> operator /(FFTWComplex<R> a, double b) { 00994 a /= b; 00995 return a; 00996 } 00997 00998 using VIGRA_CSTD::abs; 00999 01000 /// absolute value (= magnitude) 01001 template <class R> 01002 inline typename FFTWComplex<R>::NormType abs(const FFTWComplex<R> &a) 01003 { 01004 return a.magnitude(); 01005 } 01006 01007 /// pahse 01008 template <class R> 01009 inline R arg(const FFTWComplex<R> &a) 01010 { 01011 return a.phase(); 01012 } 01013 01014 /// real part 01015 template <class R> 01016 inline R real(const FFTWComplex<R> &a) 01017 { 01018 return a.real(); 01019 } 01020 01021 /// imaginary part 01022 template <class R> 01023 inline R imag(const FFTWComplex<R> &a) 01024 { 01025 return a.imag(); 01026 } 01027 01028 /// complex conjugate 01029 template <class R> 01030 inline FFTWComplex<R> conj(const FFTWComplex<R> &a) 01031 { 01032 return FFTWComplex<R>(a.re(), -a.im()); 01033 } 01034 01035 /// norm (= magnitude) 01036 template <class R> 01037 inline typename FFTWComplex<R>::NormType norm(const FFTWComplex<R> &a) 01038 { 01039 return a.magnitude(); 01040 } 01041 01042 /// squared norm (= squared magnitude) 01043 template <class R> 01044 inline typename FFTWComplex<R>::SquaredNormType squaredNorm(const FFTWComplex<R> &a) 01045 { 01046 return a.squaredMagnitude(); 01047 } 01048 01049 #define VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(fct) \ 01050 template <class R> \ 01051 inline FFTWComplex<R> fct(const FFTWComplex<R> &a) \ 01052 { \ 01053 return std::fct(reinterpret_cast<std::complex<R> const &>(a)); \ 01054 } 01055 01056 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cos) 01057 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cosh) 01058 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(exp) 01059 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(log) 01060 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(log10) 01061 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sin) 01062 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sinh) 01063 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sqrt) 01064 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tan) 01065 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tanh) 01066 01067 #undef VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION 01068 01069 template <class R> 01070 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, int e) 01071 { 01072 return std::pow(reinterpret_cast<std::complex<R> const &>(a), e); 01073 } 01074 01075 template <class R> 01076 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, R const & e) 01077 { 01078 return std::pow(reinterpret_cast<std::complex<R> const &>(a), e); 01079 } 01080 01081 template <class R> 01082 inline FFTWComplex<R> pow(const FFTWComplex<R> &a, const FFTWComplex<R> & e) 01083 { 01084 return std::pow(reinterpret_cast<std::complex<R> const &>(a), 01085 reinterpret_cast<std::complex<R> const &>(e)); 01086 } 01087 01088 template <class R> 01089 inline FFTWComplex<R> pow(R const & a, const FFTWComplex<R> &e) 01090 { 01091 return std::pow(a, reinterpret_cast<std::complex<R> const &>(e)); 01092 } 01093 01094 //@} 01095 01096 } // namespace vigra 01097 01098 namespace std { 01099 01100 template <class Real> 01101 ostream & operator<<(ostream & s, vigra::FFTWComplex<Real> const & v) 01102 { 01103 s << std::complex<Real>(v.re(), v.im()); 01104 return s; 01105 } 01106 01107 } // namespace std 01108 01109 namespace vigra { 01110 01111 /** \addtogroup StandardImageTypes 01112 */ 01113 //@{ 01114 01115 /********************************************************/ 01116 /* */ 01117 /* FFTWRealImage */ 01118 /* */ 01119 /********************************************************/ 01120 01121 /** Float (<tt>fftw_real</tt>) image. 01122 01123 The type <tt>fftw_real</tt> is defined as <tt>double</tt> (in FFTW 2 it used to be 01124 either <tt>float</tt> or <tt>double</tt>, as specified during compilation of FFTW). 01125 FFTWRealImage uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and 01126 their const counterparts to access the data. 01127 01128 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01129 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01130 Namespace: vigra 01131 */ 01132 typedef BasicImage<fftw_real> FFTWRealImage; 01133 01134 /********************************************************/ 01135 /* */ 01136 /* FFTWComplexImage */ 01137 /* */ 01138 /********************************************************/ 01139 01140 template<class R> 01141 struct IteratorTraits< 01142 BasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> > 01143 { 01144 typedef FFTWComplex<R> Type; 01145 typedef BasicImageIterator<Type, Type **> Iterator; 01146 typedef Iterator iterator; 01147 typedef BasicImageIterator<Type, Type **> mutable_iterator; 01148 typedef ConstBasicImageIterator<Type, Type **> const_iterator; 01149 typedef typename iterator::iterator_category iterator_category; 01150 typedef typename iterator::value_type value_type; 01151 typedef typename iterator::reference reference; 01152 typedef typename iterator::index_reference index_reference; 01153 typedef typename iterator::pointer pointer; 01154 typedef typename iterator::difference_type difference_type; 01155 typedef typename iterator::row_iterator row_iterator; 01156 typedef typename iterator::column_iterator column_iterator; 01157 typedef VectorAccessor<Type> default_accessor; 01158 typedef VectorAccessor<Type> DefaultAccessor; 01159 typedef VigraTrueType hasConstantStrides; 01160 }; 01161 01162 template<class R> 01163 struct IteratorTraits< 01164 ConstBasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> > 01165 { 01166 typedef FFTWComplex<R> Type; 01167 typedef ConstBasicImageIterator<Type, Type **> Iterator; 01168 typedef Iterator iterator; 01169 typedef BasicImageIterator<Type, Type **> mutable_iterator; 01170 typedef ConstBasicImageIterator<Type, Type **> const_iterator; 01171 typedef typename iterator::iterator_category iterator_category; 01172 typedef typename iterator::value_type value_type; 01173 typedef typename iterator::reference reference; 01174 typedef typename iterator::index_reference index_reference; 01175 typedef typename iterator::pointer pointer; 01176 typedef typename iterator::difference_type difference_type; 01177 typedef typename iterator::row_iterator row_iterator; 01178 typedef typename iterator::column_iterator column_iterator; 01179 typedef VectorAccessor<Type> default_accessor; 01180 typedef VectorAccessor<Type> DefaultAccessor; 01181 typedef VigraTrueType hasConstantStrides; 01182 }; 01183 01184 /** Complex (FFTWComplex) image. 01185 It uses \ref vigra::BasicImageIterator and \ref vigra::StandardAccessor and 01186 their const counterparts to access the data. 01187 01188 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01189 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01190 Namespace: vigra 01191 */ 01192 typedef BasicImage<FFTWComplex<> > FFTWComplexImage; 01193 01194 //@} 01195 01196 /********************************************************/ 01197 /* */ 01198 /* FFTWComplex-Accessors */ 01199 /* */ 01200 /********************************************************/ 01201 01202 /** \addtogroup DataAccessors 01203 */ 01204 //@{ 01205 /** \defgroup FFTWComplexAccessors Accessors for FFTWComplex 01206 01207 Encapsulate access to pixels of type FFTWComplex 01208 */ 01209 //@{ 01210 /** Encapsulate access to the the real part of a complex number. 01211 01212 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01213 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01214 Namespace: vigra 01215 */ 01216 template <class Real = double> 01217 class FFTWRealAccessor 01218 { 01219 public: 01220 01221 /// The accessor's value type. 01222 typedef Real value_type; 01223 01224 /// Read real part at iterator position. 01225 template <class ITERATOR> 01226 value_type operator()(ITERATOR const & i) const { 01227 return (*i).re(); 01228 } 01229 01230 /// Read real part at offset from iterator position. 01231 template <class ITERATOR, class DIFFERENCE> 01232 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 01233 return i[d].re(); 01234 } 01235 01236 /// Write real part at iterator position from a scalar. 01237 template <class ITERATOR> 01238 void set(value_type const & v, ITERATOR const & i) const { 01239 (*i).re()= v; 01240 } 01241 01242 /// Write real part at offset from iterator position from a scalar. 01243 template <class ITERATOR, class DIFFERENCE> 01244 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 01245 i[d].re()= v; 01246 } 01247 01248 /// Write real part at iterator position into a scalar. 01249 template <class R, class ITERATOR> 01250 void set(FFTWComplex<R> const & v, ITERATOR const & i) const { 01251 *i = v.re(); 01252 } 01253 01254 /// Write real part at offset from iterator position into a scalar. 01255 template <class R, class ITERATOR, class DIFFERENCE> 01256 void set(FFTWComplex<R> const & v, ITERATOR const & i, DIFFERENCE d) const { 01257 i[d] = v.re(); 01258 } 01259 }; 01260 01261 /** Encapsulate access to the the imaginary part of a complex number. 01262 01263 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01264 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01265 Namespace: vigra 01266 */ 01267 template <class Real = double> 01268 class FFTWImaginaryAccessor 01269 { 01270 public: 01271 /// The accessor's value type. 01272 typedef Real value_type; 01273 01274 /// Read imaginary part at iterator position. 01275 template <class ITERATOR> 01276 value_type operator()(ITERATOR const & i) const { 01277 return (*i).im(); 01278 } 01279 01280 /// Read imaginary part at offset from iterator position. 01281 template <class ITERATOR, class DIFFERENCE> 01282 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 01283 return i[d].im(); 01284 } 01285 01286 /// Write imaginary part at iterator position from a scalar. 01287 template <class ITERATOR> 01288 void set(value_type const & v, ITERATOR const & i) const { 01289 (*i).im()= v; 01290 } 01291 01292 /// Write imaginary part at offset from iterator position from a scalar. 01293 template <class ITERATOR, class DIFFERENCE> 01294 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 01295 i[d].im()= v; 01296 } 01297 01298 /// Write imaginary part at iterator position into a scalar. 01299 template <class R, class ITERATOR> 01300 void set(FFTWComplex<R> const & v, ITERATOR const & i) const { 01301 *i = v.im(); 01302 } 01303 01304 /// Write imaginary part at offset from iterator position into a scalar. 01305 template <class R, class ITERATOR, class DIFFERENCE> 01306 void set(FFTWComplex<R> const & v, ITERATOR const & i, DIFFERENCE d) const { 01307 i[d] = v.im(); 01308 } 01309 }; 01310 01311 /** Write a real number into a complex one. The imaginary part is set 01312 to 0. 01313 01314 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01315 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01316 Namespace: vigra 01317 */ 01318 template <class Real = double> 01319 class FFTWWriteRealAccessor 01320 : public FFTWRealAccessor<Real> 01321 { 01322 public: 01323 /// The accessor's value type. 01324 typedef Real value_type; 01325 01326 /** Write real number at iterator position. Set imaginary part 01327 to 0. 01328 */ 01329 template <class ITERATOR> 01330 void set(value_type const & v, ITERATOR const & i) const { 01331 (*i).re()= v; 01332 (*i).im()= 0; 01333 } 01334 01335 /** Write real number at offset from iterator position. Set imaginary part 01336 to 0. 01337 */ 01338 template <class ITERATOR, class DIFFERENCE> 01339 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 01340 i[d].re()= v; 01341 i[d].im()= 0; 01342 } 01343 }; 01344 01345 /** Calculate squared magnitude of complex number on the fly. 01346 01347 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01348 Namespace: vigra 01349 */ 01350 template <class Real = double> 01351 class FFTWSquaredMagnitudeAccessor 01352 { 01353 public: 01354 /// The accessor's value type. 01355 typedef Real value_type; 01356 01357 /// Read squared magnitude at iterator position. 01358 template <class ITERATOR> 01359 value_type operator()(ITERATOR const & i) const { 01360 return (*i).squaredMagnitude(); 01361 } 01362 01363 /// Read squared magnitude at offset from iterator position. 01364 template <class ITERATOR, class DIFFERENCE> 01365 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 01366 return (i[d]).squaredMagnitude(); 01367 } 01368 }; 01369 01370 /** Calculate magnitude of complex number on the fly. 01371 01372 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01373 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01374 Namespace: vigra 01375 */ 01376 template <class Real = double> 01377 class FFTWMagnitudeAccessor 01378 { 01379 public: 01380 /// The accessor's value type. 01381 typedef Real value_type; 01382 01383 /// Read magnitude at iterator position. 01384 template <class ITERATOR> 01385 value_type operator()(ITERATOR const & i) const { 01386 return (*i).magnitude(); 01387 } 01388 01389 /// Read magnitude at offset from iterator position. 01390 template <class ITERATOR, class DIFFERENCE> 01391 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 01392 return (i[d]).magnitude(); 01393 } 01394 }; 01395 01396 /** Calculate phase of complex number on the fly. 01397 01398 <b>\#include</b> <vigra/fftw3.hxx> (for FFTW 3) or<br> 01399 <b>\#include</b> <vigra/fftw.hxx> (for deprecated FFTW 2)<br> 01400 Namespace: vigra 01401 */ 01402 template <class Real = double> 01403 class FFTWPhaseAccessor 01404 { 01405 public: 01406 /// The accessor's value type. 01407 typedef Real value_type; 01408 01409 /// Read phase at iterator position. 01410 template <class ITERATOR> 01411 value_type operator()(ITERATOR const & i) const { 01412 return (*i).phase(); 01413 } 01414 01415 /// Read phase at offset from iterator position. 01416 template <class ITERATOR, class DIFFERENCE> 01417 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 01418 return (i[d]).phase(); 01419 } 01420 }; 01421 01422 //@} 01423 //@} 01424 01425 /********************************************************/ 01426 /* */ 01427 /* Fourier Transform */ 01428 /* */ 01429 /********************************************************/ 01430 01431 /* \addtogroup FourierTransform Fast Fourier Transform 01432 01433 This documentation describes the VIGRA interface to FFTW version 3. The interface 01434 to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated. 01435 01436 VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier 01437 Transform</a> package to perform Fourier transformations. VIGRA 01438 provides a wrapper for FFTW's complex number type (FFTWComplex), 01439 but FFTW's functions are used verbatim. If the image is stored as 01440 a FFTWComplexImage, the simplest call to an FFT function is like this: 01441 01442 \code 01443 vigra::FFTWComplexImage spatial(width,height), fourier(width,height); 01444 ... // fill image with data 01445 01446 // create a plan with estimated performance optimization 01447 fftw_plan forwardPlan = fftw_plan_dft_2d(height, width, 01448 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), 01449 FFTW_FORWARD, FFTW_ESTIMATE ); 01450 // calculate FFT (this can be repeated as often as needed, 01451 // with fresh data written into the source array) 01452 fftw_execute(forwardPlan); 01453 01454 // release the plan memory 01455 fftw_destroy_plan(forwardPlan); 01456 01457 // likewise for the inverse transform 01458 fftw_plan backwardPlan = fftw_plan_dft_2d(height, width, 01459 (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(), 01460 FFTW_BACKWARD, FFTW_ESTIMATE); 01461 fftw_execute(backwardPlan); 01462 fftw_destroy_plan(backwardPlan); 01463 01464 // do not forget to normalize the result according to the image size 01465 transformImage(srcImageRange(spatial), destImage(spatial), 01466 std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height)); 01467 \endcode 01468 01469 Note that in the creation of a plan, the height must be given 01470 first. Note also that <TT>spatial.begin()</TT> may only be passed 01471 to <TT>fftw_plan_dft_2d</TT> if the transform shall be applied to the 01472 entire image. When you want to restrict operation to an ROI, you 01473 can create a copy of the ROI in an image of appropriate size, or 01474 you may use the Guru interface to FFTW. 01475 01476 More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>. 01477 01478 FFTW produces fourier images that have the DC component (the 01479 origin of the Fourier space) in the upper left corner. Often, one 01480 wants the origin in the center of the image, so that frequencies 01481 always increase towards the border of the image. This can be 01482 achieved by calling \ref moveDCToCenter(). The inverse 01483 transformation is done by \ref moveDCToUpperLeft(). 01484 01485 <b>\#include</b> <vigra/fftw3.hxx><br> 01486 Namespace: vigra 01487 */ 01488 01489 /** \addtogroup FourierTransform Fast Fourier Transform 01490 01491 VIGRA provides a powerful C++ API for the popular <a href="http://www.fftw.org/">FFTW library</a> 01492 for fast Fourier transforms. There are two versions of the API: an older one based on image 01493 iterators (and therefore restricted to 2D) and a new one based on \ref MultiArrayView that 01494 works for arbitrary dimensions. In addition, the functions \ref convolveFFT() and 01495 \ref applyFourierFilter() provide an easy-to-use interface for FFT-based convolution, 01496 a major application of Fourier transforms. 01497 */ 01498 //@{ 01499 01500 /********************************************************/ 01501 /* */ 01502 /* moveDCToCenter */ 01503 /* */ 01504 /********************************************************/ 01505 01506 /** \brief Rearrange the quadrants of a Fourier image so that the origin is 01507 in the image center. 01508 01509 FFTW produces fourier images where the DC component (origin of 01510 fourier space) is located in the upper left corner of the 01511 image. The quadrants are placed like this (using a 4x4 image for 01512 example): 01513 01514 \code 01515 DC 4 3 3 01516 4 4 3 3 01517 1 1 2 2 01518 1 1 2 2 01519 \endcode 01520 01521 After applying the function, the quadrants are at their usual places: 01522 01523 \code 01524 2 2 1 1 01525 2 2 1 1 01526 3 3 DC 4 01527 3 3 4 4 01528 \endcode 01529 01530 This transformation can be reversed by \ref moveDCToUpperLeft(). 01531 Note that the 2D versions of this transformation must not be executed in place - input 01532 and output images must be different. In contrast, the nD version (with MultiArrayView 01533 argument) always works in-place. 01534 01535 <b> Declarations:</b> 01536 01537 use MultiArrayView (this works in-place, with arbitrary dimension N): 01538 \code 01539 namespace vigra { 01540 template <unsigned int N, class T, class Stride> 01541 void moveDCToCenter(MultiArrayView<N, T, Stride> a); 01542 } 01543 \endcode 01544 01545 pass iterators explicitly (2D only, not in-place): 01546 \code 01547 namespace vigra { 01548 template <class SrcImageIterator, class SrcAccessor, 01549 class DestImageIterator, class DestAccessor> 01550 void moveDCToCenter(SrcImageIterator src_upperleft, 01551 SrcImageIterator src_lowerright, SrcAccessor sa, 01552 DestImageIterator dest_upperleft, DestAccessor da); 01553 } 01554 \endcode 01555 01556 01557 use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place): 01558 \code 01559 namespace vigra { 01560 template <class SrcImageIterator, class SrcAccessor, 01561 class DestImageIterator, class DestAccessor> 01562 void moveDCToCenter( 01563 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01564 pair<DestImageIterator, DestAccessor> dest); 01565 } 01566 \endcode 01567 01568 <b> Usage:</b> 01569 01570 <b>\#include</b> <vigra/fftw3.hxx> (for 2D variants) <br> 01571 <b>\#include</b> <vigra/multi_fft.hxx> (for nD variants) <br> 01572 Namespace: vigra 01573 01574 \code 01575 vigra::FFTWComplexImage spatial(width,height), fourier(width,height); 01576 ... // fill image with data 01577 01578 // create a plan with estimated performance optimization 01579 fftw_plan forwardPlan = fftw_plan_dft_2d(height, width, 01580 (fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(), 01581 FFTW_FORWARD, FFTW_ESTIMATE ); 01582 // calculate FFT 01583 fftw_execute(forwardPlan); 01584 01585 vigra::FFTWComplexImage rearrangedFourier(width, height); 01586 moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier)); 01587 01588 // delete the plan 01589 fftw_destroy_plan(forwardPlan); 01590 \endcode 01591 */ 01592 doxygen_overloaded_function(template <...> void moveDCToCenter) 01593 01594 template <class SrcImageIterator, class SrcAccessor, 01595 class DestImageIterator, class DestAccessor> 01596 void moveDCToCenter(SrcImageIterator src_upperleft, 01597 SrcImageIterator src_lowerright, SrcAccessor sa, 01598 DestImageIterator dest_upperleft, DestAccessor da) 01599 { 01600 int w = int(src_lowerright.x - src_upperleft.x); 01601 int h = int(src_lowerright.y - src_upperleft.y); 01602 int w1 = w/2; 01603 int h1 = h/2; 01604 int w2 = (w+1)/2; 01605 int h2 = (h+1)/2; 01606 01607 // 2. Quadrant zum 4. 01608 copyImage(srcIterRange(src_upperleft, 01609 src_upperleft + Diff2D(w2, h2), sa), 01610 destIter (dest_upperleft + Diff2D(w1, h1), da)); 01611 01612 // 4. Quadrant zum 2. 01613 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2), 01614 src_lowerright, sa), 01615 destIter (dest_upperleft, da)); 01616 01617 // 1. Quadrant zum 3. 01618 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0), 01619 src_upperleft + Diff2D(w, h2), sa), 01620 destIter (dest_upperleft + Diff2D(0, h1), da)); 01621 01622 // 3. Quadrant zum 1. 01623 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2), 01624 src_upperleft + Diff2D(w2, h), sa), 01625 destIter (dest_upperleft + Diff2D(w1, 0), da)); 01626 } 01627 01628 template <class SrcImageIterator, class SrcAccessor, 01629 class DestImageIterator, class DestAccessor> 01630 inline void moveDCToCenter( 01631 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01632 pair<DestImageIterator, DestAccessor> dest) 01633 { 01634 moveDCToCenter(src.first, src.second, src.third, 01635 dest.first, dest.second); 01636 } 01637 01638 /********************************************************/ 01639 /* */ 01640 /* moveDCToUpperLeft */ 01641 /* */ 01642 /********************************************************/ 01643 01644 /** \brief Rearrange the quadrants of a Fourier image so that the origin is 01645 in the image's upper left. 01646 01647 This function is the inversion of \ref moveDCToCenter(). See there 01648 for a detailed description and usage examples. 01649 01650 <b> Declarations:</b> 01651 01652 use MultiArrayView (this works in-place, with arbitrary dimension N): 01653 \code 01654 namespace vigra { 01655 template <unsigned int N, class T, class Stride> 01656 void moveDCToUpperLeft(MultiArrayView<N, T, Stride> a); 01657 } 01658 \endcode 01659 01660 pass iterators explicitly (2D only, not in-place): 01661 \code 01662 namespace vigra { 01663 template <class SrcImageIterator, class SrcAccessor, 01664 class DestImageIterator, class DestAccessor> 01665 void moveDCToUpperLeft(SrcImageIterator src_upperleft, 01666 SrcImageIterator src_lowerright, SrcAccessor sa, 01667 DestImageIterator dest_upperleft, DestAccessor da); 01668 } 01669 \endcode 01670 01671 01672 use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, not in-place): 01673 \code 01674 namespace vigra { 01675 template <class SrcImageIterator, class SrcAccessor, 01676 class DestImageIterator, class DestAccessor> 01677 void moveDCToUpperLeft( 01678 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01679 pair<DestImageIterator, DestAccessor> dest); 01680 } 01681 \endcode 01682 */ 01683 doxygen_overloaded_function(template <...> void moveDCToUpperLeft) 01684 01685 template <class SrcImageIterator, class SrcAccessor, 01686 class DestImageIterator, class DestAccessor> 01687 void moveDCToUpperLeft(SrcImageIterator src_upperleft, 01688 SrcImageIterator src_lowerright, SrcAccessor sa, 01689 DestImageIterator dest_upperleft, DestAccessor da) 01690 { 01691 int w = int(src_lowerright.x - src_upperleft.x); 01692 int h = int(src_lowerright.y - src_upperleft.y); 01693 int w2 = w/2; 01694 int h2 = h/2; 01695 int w1 = (w+1)/2; 01696 int h1 = (h+1)/2; 01697 01698 // 2. Quadrant zum 4. 01699 copyImage(srcIterRange(src_upperleft, 01700 src_upperleft + Diff2D(w2, h2), sa), 01701 destIter (dest_upperleft + Diff2D(w1, h1), da)); 01702 01703 // 4. Quadrant zum 2. 01704 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2), 01705 src_lowerright, sa), 01706 destIter (dest_upperleft, da)); 01707 01708 // 1. Quadrant zum 3. 01709 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0), 01710 src_upperleft + Diff2D(w, h2), sa), 01711 destIter (dest_upperleft + Diff2D(0, h1), da)); 01712 01713 // 3. Quadrant zum 1. 01714 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2), 01715 src_upperleft + Diff2D(w2, h), sa), 01716 destIter (dest_upperleft + Diff2D(w1, 0), da)); 01717 } 01718 01719 template <class SrcImageIterator, class SrcAccessor, 01720 class DestImageIterator, class DestAccessor> 01721 inline void moveDCToUpperLeft( 01722 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01723 pair<DestImageIterator, DestAccessor> dest) 01724 { 01725 moveDCToUpperLeft(src.first, src.second, src.third, 01726 dest.first, dest.second); 01727 } 01728 01729 namespace detail { 01730 01731 template <class T> 01732 void 01733 fourierTransformImpl(FFTWComplexImage::const_traverser sul, 01734 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01735 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest, T sign) 01736 { 01737 int w = int(slr.x - sul.x); 01738 int h = int(slr.y - sul.y); 01739 01740 FFTWComplexImage sworkImage, dworkImage; 01741 01742 fftw_complex * srcPtr = (fftw_complex *)(&*sul); 01743 fftw_complex * destPtr = (fftw_complex *)(&*dul); 01744 01745 // test for right memory layout (fftw expects a 2*width*height floats array) 01746 if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1)))) 01747 { 01748 sworkImage.resize(w, h); 01749 copyImage(srcIterRange(sul, slr, src), destImage(sworkImage)); 01750 srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft())); 01751 } 01752 if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1)))) 01753 { 01754 dworkImage.resize(w, h); 01755 destPtr = (fftw_complex *)(&(*dworkImage.upperLeft())); 01756 } 01757 01758 fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE ); 01759 fftw_execute(plan); 01760 fftw_destroy_plan(plan); 01761 01762 if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1)))) 01763 { 01764 copyImage(srcImageRange(dworkImage), destIter(dul, dest)); 01765 } 01766 } 01767 01768 } // namespace detail 01769 01770 /********************************************************/ 01771 /* */ 01772 /* fourierTransform */ 01773 /* */ 01774 /********************************************************/ 01775 01776 /** \brief Compute forward and inverse Fourier transforms. 01777 01778 The array referring to the spatial domain (i.e. the input in a forward transform, 01779 and the output in an inverse transform) may be scalar or complex. The array representing 01780 the frequency domain (i.e. output for forward transform, input for inverse transform) 01781 must always be complex. 01782 01783 The new implementations (those using MultiArrayView arguments) perform a normalized transform, 01784 whereas the old ones (using 2D iterators or argument objects) perform an un-normalized 01785 transform (i.e. the result of the inverse transform is scaled by the number of pixels). 01786 01787 In general, input and output arrays must have the same shape, with the exception of the 01788 special <a href="http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html">R2C 01789 and C2R modes</a> defined by FFTW. 01790 01791 The R2C transform reduces the redundancy in the Fourier representation of a real-valued signal: 01792 Since the Fourier representation of a real signal is symmetric, about half of the Fourier coefficients 01793 can simply be dropped. By convention, this reduction is applied to the first (innermost) dimension, 01794 such that <tt>fourier.shape(0) == spatial.shape(0)/2 + 1</tt> holds. The correct frequency domain 01795 shape can be conveniently computed by means of the function \ref fftwCorrespondingShapeR2C(). 01796 01797 Note that your program must always link against <tt>libfftw3</tt>. If you want to compute Fourier 01798 transforms for <tt>float</tt> or <tt>long double</tt> arrays, you must <i>additionally</i> link against <tt>libfftw3f</tt> and <tt>libfftw3l</tt> respectively. (Old-style functions only support <tt>double</tt>). 01799 01800 The Fourier transform functions internally create <a href="http://www.fftw.org/doc/Using-Plans.html">FFTW plans</a> 01801 which control the algorithm details. The plans are creates with the flag <tt>FFTW_ESTIMATE</tt>, i.e. 01802 optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning, 01803 you can use the class \ref FFTWPlan. 01804 01805 <b> Declarations:</b> 01806 01807 use complex-valued MultiArrayView arguments (this works for arbitrary dimension N): 01808 \code 01809 namespace vigra { 01810 template <unsigned int N, class Real, class C1, class C2> 01811 void 01812 fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01813 MultiArrayView<N, FFTWComplex<Real>, C2> out); 01814 01815 template <unsigned int N, class Real, class C1, class C2> 01816 void 01817 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01818 MultiArrayView<N, FFTWComplex<Real>, C2> out); 01819 } 01820 \endcode 01821 01822 use real-valued MultiArrayView in the spatial domain, complex-valued MultiArrayView 01823 in the frequency domain (this works for arbitrary dimension N, and also supports 01824 the R2C and C2R transform, depending on the array shape in the frequency domain): 01825 \code 01826 namespace vigra { 01827 template <unsigned int N, class Real, class C1, class C2> 01828 void 01829 fourierTransform(MultiArrayView<N, Real, C1> in, 01830 MultiArrayView<N, FFTWComplex<Real>, C2> out); 01831 01832 template <unsigned int N, class Real, class C1, class C2> 01833 void 01834 fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in, 01835 MultiArrayView<N, Real, C2> out); 01836 } 01837 \endcode 01838 01839 pass iterators explicitly (2D only, double only): 01840 \code 01841 namespace vigra { 01842 template <class SrcImageIterator, class SrcAccessor> 01843 void fourierTransform(SrcImageIterator srcUpperLeft, 01844 SrcImageIterator srcLowerRight, SrcAccessor src, 01845 FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest); 01846 01847 void 01848 fourierTransformInverse(FFTWComplexImage::const_traverser sul, 01849 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01850 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest) 01851 } 01852 \endcode 01853 01854 use argument objects in conjunction with \ref ArgumentObjectFactories (2D only, double only): 01855 \code 01856 namespace vigra { 01857 template <class SrcImageIterator, class SrcAccessor> 01858 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01859 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest); 01860 01861 void 01862 fourierTransformInverse(triple<FFTWComplexImage::const_traverser, 01863 FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src, 01864 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest); 01865 } 01866 \endcode 01867 01868 <b> Usage:</b> 01869 01870 <b>\#include</b> <vigra/fftw3.hxx> (old-style 2D variants)<br> 01871 <b>\#include</b> <vigra/multi_fft.hxx> (new-style nD variants)<br> 01872 Namespace: vigra 01873 01874 old-style example (using FFTWComplexImage): 01875 \code 01876 // compute complex Fourier transform of a real image, old-style implementation 01877 vigra::DImage src(w, h); 01878 vigra::FFTWComplexImage fourier(w, h); 01879 01880 fourierTransform(srcImageRange(src), destImage(fourier)); 01881 01882 // compute inverse Fourier transform 01883 // note that both source and destination image must be of type vigra::FFTWComplexImage 01884 vigra::FFTWComplexImage inverseFourier(w, h); 01885 01886 fourierTransformInverse(srcImageRange(fourier), destImage(inverseFourier)); 01887 \endcode 01888 01889 new-style examples (using MultiArray): 01890 \code 01891 // compute Fourier transform of a real array, using the R2C algorithm 01892 MultiArray<2, double> src(Shape2(w, h)); 01893 MultiArray<2, FFTWComplex<double> > fourier(fftwCorrespondingShapeR2C(src.shape())); 01894 01895 fourierTransform(src, fourier); 01896 01897 // compute inverse Fourier transform, using the C2R algorithm 01898 MultiArray<2, double> dest(src.shape()); 01899 fourierTransformInverse(fourier, dest); 01900 \endcode 01901 01902 \code 01903 // compute Fourier transform of a real array with standard algorithm 01904 MultiArray<2, double> src(Shape2(w, h)); 01905 MultiArray<2, FFTWComplex<double> > fourier(src.shape()); 01906 01907 fourierTransform(src, fourier); 01908 01909 // compute inverse Fourier transform, using the C2R algorithm 01910 MultiArray<2, double> dest(src.shape()); 01911 fourierTransformInverse(fourier, dest); 01912 \endcode 01913 Complex input arrays are handled in the same way. 01914 01915 */ 01916 doxygen_overloaded_function(template <...> void fourierTransform) 01917 01918 inline void 01919 fourierTransform(FFTWComplexImage::const_traverser sul, 01920 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01921 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest) 01922 { 01923 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD); 01924 } 01925 01926 template <class SrcImageIterator, class SrcAccessor> 01927 void fourierTransform(SrcImageIterator srcUpperLeft, 01928 SrcImageIterator srcLowerRight, SrcAccessor sa, 01929 FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor da) 01930 { 01931 // copy real input images into a complex one... 01932 int w= srcLowerRight.x - srcUpperLeft.x; 01933 int h= srcLowerRight.y - srcUpperLeft.y; 01934 01935 FFTWComplexImage workImage(w, h); 01936 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01937 destImage(workImage, FFTWWriteRealAccessor<>())); 01938 01939 // ...and call the complex -> complex version of the algorithm 01940 FFTWComplexImage const & cworkImage = workImage; 01941 fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01942 destUpperLeft, da); 01943 } 01944 01945 template <class SrcImageIterator, class SrcAccessor> 01946 inline 01947 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01948 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest) 01949 { 01950 fourierTransform(src.first, src.second, src.third, dest.first, dest.second); 01951 } 01952 01953 /** \brief Compute inverse Fourier transforms. 01954 01955 See \ref fourierTransform() for details. 01956 */ 01957 doxygen_overloaded_function(template <...> void fourierTransformInverse) 01958 01959 inline void 01960 fourierTransformInverse(FFTWComplexImage::const_traverser sul, 01961 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01962 FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest) 01963 { 01964 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD); 01965 } 01966 01967 template <class DestImageIterator, class DestAccessor> 01968 void fourierTransformInverse(FFTWComplexImage::const_traverser sul, 01969 FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, 01970 DestImageIterator dul, DestAccessor dest) 01971 { 01972 int w = slr.x - sul.x; 01973 int h = slr.y - sul.y; 01974 01975 FFTWComplexImage workImage(w, h); 01976 fourierTransformInverse(sul, slr, src, workImage.upperLeft(), workImage.accessor()); 01977 copyImage(srcImageRange(workImage), destIter(dul, dest)); 01978 } 01979 01980 01981 template <class DestImageIterator, class DestAccessor> 01982 inline void 01983 fourierTransformInverse(triple<FFTWComplexImage::const_traverser, 01984 FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src, 01985 pair<DestImageIterator, DestAccessor> dest) 01986 { 01987 fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second); 01988 } 01989 01990 01991 01992 /********************************************************/ 01993 /* */ 01994 /* applyFourierFilter */ 01995 /* */ 01996 /********************************************************/ 01997 01998 /** \brief Apply a filter (defined in the frequency domain) to an image. 01999 02000 After transferring the image into the frequency domain, it is 02001 multiplied pixel-wise with the filter and transformed back. The 02002 result is put into the given destination image which must have the right size. 02003 The result will be normalized to compensate for the two FFTs. 02004 02005 If the destination image is scalar, only the real part of the result image is 02006 retained. In this case, you are responsible for choosing a filter image 02007 which ensures a zero imaginary part of the result (e.g. use a real, even symmetric 02008 filter image, or a purely imaginary, odd symmetric one). 02009 02010 The DC entry of the filter must be in the upper left, which is the 02011 position where FFTW expects it (see \ref moveDCToUpperLeft()). 02012 02013 See also \ref convolveFFT() for corresponding functionality on the basis of the 02014 \ref MultiArrayView interface. 02015 02016 <b> Declarations:</b> 02017 02018 pass arguments explicitly: 02019 \code 02020 namespace vigra { 02021 template <class SrcImageIterator, class SrcAccessor, 02022 class FilterImageIterator, class FilterAccessor, 02023 class DestImageIterator, class DestAccessor> 02024 void applyFourierFilter(SrcImageIterator srcUpperLeft, 02025 SrcImageIterator srcLowerRight, SrcAccessor sa, 02026 FilterImageIterator filterUpperLeft, FilterAccessor fa, 02027 DestImageIterator destUpperLeft, DestAccessor da); 02028 } 02029 \endcode 02030 02031 use argument objects in conjunction with \ref ArgumentObjectFactories : 02032 \code 02033 namespace vigra { 02034 template <class SrcImageIterator, class SrcAccessor, 02035 class FilterImageIterator, class FilterAccessor, 02036 class DestImageIterator, class DestAccessor> 02037 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 02038 pair<FilterImageIterator, FilterAccessor> filter, 02039 pair<DestImageIterator, DestAccessor> dest); 02040 } 02041 \endcode 02042 02043 <b> Usage:</b> 02044 02045 <b>\#include</b> <vigra/fftw3.hxx><br> 02046 Namespace: vigra 02047 02048 \code 02049 // create a Gaussian filter in Fourier space 02050 vigra::FImage gaussFilter(w, h), filter(w, h); 02051 for(int y=0; y<h; ++y) 02052 for(int x=0; x<w; ++x) 02053 { 02054 xx = float(x - w / 2) / w; 02055 yy = float(y - h / 2) / h; 02056 02057 gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale); 02058 } 02059 02060 // applyFourierFilter() expects the filter's DC in the upper left 02061 moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter)); 02062 02063 vigra::FFTWComplexImage result(w, h); 02064 02065 vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result); 02066 \endcode 02067 02068 For inspection of the result, \ref FFTWMagnitudeAccessor might be 02069 useful. If you want to apply the same filter repeatedly, it may be more 02070 efficient to use the FFTW functions directly with FFTW plans optimized 02071 for good performance. 02072 */ 02073 doxygen_overloaded_function(template <...> void applyFourierFilter) 02074 02075 template <class SrcImageIterator, class SrcAccessor, 02076 class FilterImageIterator, class FilterAccessor, 02077 class DestImageIterator, class DestAccessor> 02078 void applyFourierFilter(SrcImageIterator srcUpperLeft, 02079 SrcImageIterator srcLowerRight, SrcAccessor sa, 02080 FilterImageIterator filterUpperLeft, FilterAccessor fa, 02081 DestImageIterator destUpperLeft, DestAccessor da) 02082 { 02083 // copy real input images into a complex one... 02084 int w = int(srcLowerRight.x - srcUpperLeft.x); 02085 int h = int(srcLowerRight.y - srcUpperLeft.y); 02086 02087 FFTWComplexImage workImage(w, h); 02088 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 02089 destImage(workImage, FFTWWriteRealAccessor<>())); 02090 02091 // ...and call the impl 02092 FFTWComplexImage const & cworkImage = workImage; 02093 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 02094 filterUpperLeft, fa, 02095 destUpperLeft, da); 02096 } 02097 02098 template <class FilterImageIterator, class FilterAccessor, 02099 class DestImageIterator, class DestAccessor> 02100 inline 02101 void applyFourierFilter( 02102 FFTWComplexImage::const_traverser srcUpperLeft, 02103 FFTWComplexImage::const_traverser srcLowerRight, 02104 FFTWComplexImage::ConstAccessor sa, 02105 FilterImageIterator filterUpperLeft, FilterAccessor fa, 02106 DestImageIterator destUpperLeft, DestAccessor da) 02107 { 02108 int w = srcLowerRight.x - srcUpperLeft.x; 02109 int h = srcLowerRight.y - srcUpperLeft.y; 02110 02111 // test for right memory layout (fftw expects a 2*width*height floats array) 02112 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 02113 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa, 02114 filterUpperLeft, fa, 02115 destUpperLeft, da); 02116 else 02117 { 02118 FFTWComplexImage workImage(w, h); 02119 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 02120 destImage(workImage)); 02121 02122 FFTWComplexImage const & cworkImage = workImage; 02123 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 02124 filterUpperLeft, fa, 02125 destUpperLeft, da); 02126 } 02127 } 02128 02129 template <class SrcImageIterator, class SrcAccessor, 02130 class FilterImageIterator, class FilterAccessor, 02131 class DestImageIterator, class DestAccessor> 02132 inline 02133 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 02134 pair<FilterImageIterator, FilterAccessor> filter, 02135 pair<DestImageIterator, DestAccessor> dest) 02136 { 02137 applyFourierFilter(src.first, src.second, src.third, 02138 filter.first, filter.second, 02139 dest.first, dest.second); 02140 } 02141 02142 template <class FilterImageIterator, class FilterAccessor, 02143 class DestImageIterator, class DestAccessor> 02144 void applyFourierFilterImpl( 02145 FFTWComplexImage::const_traverser srcUpperLeft, 02146 FFTWComplexImage::const_traverser srcLowerRight, 02147 FFTWComplexImage::ConstAccessor, 02148 FilterImageIterator filterUpperLeft, FilterAccessor fa, 02149 DestImageIterator destUpperLeft, DestAccessor da) 02150 { 02151 int w = int(srcLowerRight.x - srcUpperLeft.x); 02152 int h = int(srcLowerRight.y - srcUpperLeft.y); 02153 02154 FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft); 02155 02156 // FFT from srcImage to complexResultImg 02157 fftw_plan forwardPlan= 02158 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft), 02159 (fftw_complex *)complexResultImg.begin(), 02160 FFTW_FORWARD, FFTW_ESTIMATE ); 02161 fftw_execute(forwardPlan); 02162 fftw_destroy_plan(forwardPlan); 02163 02164 // convolve in freq. domain (in complexResultImg) 02165 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa), 02166 destImage(complexResultImg), std::multiplies<FFTWComplex<> >()); 02167 02168 // FFT back into spatial domain (inplace in complexResultImg) 02169 fftw_plan backwardPlan= 02170 fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(), 02171 (fftw_complex *)complexResultImg.begin(), 02172 FFTW_BACKWARD, FFTW_ESTIMATE); 02173 fftw_execute(backwardPlan); 02174 fftw_destroy_plan(backwardPlan); 02175 02176 typedef typename 02177 NumericTraits<typename DestAccessor::value_type>::isScalar 02178 isScalarResult; 02179 02180 // normalization (after FFTs), maybe stripping imaginary part 02181 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da, 02182 isScalarResult()); 02183 } 02184 02185 template <class DestImageIterator, class DestAccessor> 02186 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage, 02187 DestImageIterator destUpperLeft, 02188 DestAccessor da, 02189 VigraFalseType) 02190 { 02191 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 02192 02193 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 02194 { 02195 DestImageIterator dIt= destUpperLeft; 02196 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 02197 { 02198 da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0); 02199 da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1); 02200 } 02201 } 02202 } 02203 02204 inline 02205 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 02206 FFTWComplexImage::traverser destUpperLeft, 02207 FFTWComplexImage::Accessor da, 02208 VigraFalseType) 02209 { 02210 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da), 02211 linearIntensityTransform<FFTWComplex<> >(1.0/(srcImage.width() * srcImage.height()))); 02212 } 02213 02214 template <class DestImageIterator, class DestAccessor> 02215 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 02216 DestImageIterator destUpperLeft, 02217 DestAccessor da, 02218 VigraTrueType) 02219 { 02220 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 02221 02222 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 02223 { 02224 DestImageIterator dIt= destUpperLeft; 02225 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 02226 da.set(srcImage(x, y).re()*normFactor, dIt); 02227 } 02228 } 02229 02230 /**********************************************************/ 02231 /* */ 02232 /* applyFourierFilterFamily */ 02233 /* */ 02234 /**********************************************************/ 02235 02236 /** \brief Apply an array of filters (defined in the frequency domain) to an image. 02237 02238 This provides the same functionality as \ref applyFourierFilter(), 02239 but applying several filters at once allows to avoid 02240 repeated Fourier transforms of the source image. 02241 02242 Filters and result images must be stored in \ref vigra::ImageArray data 02243 structures. In contrast to \ref applyFourierFilter(), this function adjusts 02244 the size of the result images and the the length of the array. 02245 02246 <b> Declarations:</b> 02247 02248 pass arguments explicitly: 02249 \code 02250 namespace vigra { 02251 template <class SrcImageIterator, class SrcAccessor, class FilterType> 02252 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 02253 SrcImageIterator srcLowerRight, SrcAccessor sa, 02254 const ImageArray<FilterType> &filters, 02255 ImageArray<FFTWComplexImage> &results) 02256 } 02257 \endcode 02258 02259 use argument objects in conjunction with \ref ArgumentObjectFactories : 02260 \code 02261 namespace vigra { 02262 template <class SrcImageIterator, class SrcAccessor, class FilterType> 02263 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 02264 const ImageArray<FilterType> &filters, 02265 ImageArray<FFTWComplexImage> &results) 02266 } 02267 \endcode 02268 02269 <b> Usage:</b> 02270 02271 <b>\#include</b> <vigra/fftw3.hxx><br> 02272 Namespace: vigra 02273 02274 \code 02275 // assuming the presence of a real-valued image named "image" to 02276 // be filtered in this example 02277 02278 vigra::ImageArray<vigra::FImage> filters(16, image.size()); 02279 02280 for (int i=0; i<filters.size(); i++) 02281 // create some meaningful filters here 02282 createMyFilterOfScale(i, destImage(filters[i])); 02283 02284 vigra::ImageArray<vigra::FFTWComplexImage> results(); 02285 02286 vigra::applyFourierFilterFamily(srcImageRange(image), filters, results); 02287 \endcode 02288 */ 02289 doxygen_overloaded_function(template <...> void applyFourierFilterFamily) 02290 02291 template <class SrcImageIterator, class SrcAccessor, 02292 class FilterType, class DestImage> 02293 inline 02294 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 02295 const ImageArray<FilterType> &filters, 02296 ImageArray<DestImage> &results) 02297 { 02298 applyFourierFilterFamily(src.first, src.second, src.third, 02299 filters, results); 02300 } 02301 02302 template <class SrcImageIterator, class SrcAccessor, 02303 class FilterType, class DestImage> 02304 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 02305 SrcImageIterator srcLowerRight, SrcAccessor sa, 02306 const ImageArray<FilterType> &filters, 02307 ImageArray<DestImage> &results) 02308 { 02309 int w = int(srcLowerRight.x - srcUpperLeft.x); 02310 int h = int(srcLowerRight.y - srcUpperLeft.y); 02311 02312 FFTWComplexImage workImage(w, h); 02313 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 02314 destImage(workImage, FFTWWriteRealAccessor<>())); 02315 02316 FFTWComplexImage const & cworkImage = workImage; 02317 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 02318 filters, results); 02319 } 02320 02321 template <class FilterType, class DestImage> 02322 inline 02323 void applyFourierFilterFamily( 02324 FFTWComplexImage::const_traverser srcUpperLeft, 02325 FFTWComplexImage::const_traverser srcLowerRight, 02326 FFTWComplexImage::ConstAccessor sa, 02327 const ImageArray<FilterType> &filters, 02328 ImageArray<DestImage> &results) 02329 { 02330 int w= srcLowerRight.x - srcUpperLeft.x; 02331 02332 // test for right memory layout (fftw expects a 2*width*height floats array) 02333 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 02334 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa, 02335 filters, results); 02336 else 02337 { 02338 int h = srcLowerRight.y - srcUpperLeft.y; 02339 FFTWComplexImage workImage(w, h); 02340 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 02341 destImage(workImage)); 02342 02343 FFTWComplexImage const & cworkImage = workImage; 02344 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 02345 filters, results); 02346 } 02347 } 02348 02349 template <class FilterType, class DestImage> 02350 void applyFourierFilterFamilyImpl( 02351 FFTWComplexImage::const_traverser srcUpperLeft, 02352 FFTWComplexImage::const_traverser srcLowerRight, 02353 FFTWComplexImage::ConstAccessor sa, 02354 const ImageArray<FilterType> &filters, 02355 ImageArray<DestImage> &results) 02356 { 02357 // FIXME: sa is not used 02358 // (maybe check if StandardAccessor, else copy?) 02359 02360 // make sure the filter images have the right dimensions 02361 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(), 02362 "applyFourierFilterFamily called with src image size != filters.imageSize()!"); 02363 02364 // make sure the result image array has the right dimensions 02365 results.resize(filters.size()); 02366 results.resizeImages(filters.imageSize()); 02367 02368 // FFT from srcImage to freqImage 02369 int w = int(srcLowerRight.x - srcUpperLeft.x); 02370 int h = int(srcLowerRight.y - srcUpperLeft.y); 02371 02372 FFTWComplexImage freqImage(w, h); 02373 FFTWComplexImage result(w, h); 02374 02375 fftw_plan forwardPlan= 02376 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft), 02377 (fftw_complex *)freqImage.begin(), 02378 FFTW_FORWARD, FFTW_ESTIMATE ); 02379 fftw_execute(forwardPlan); 02380 fftw_destroy_plan(forwardPlan); 02381 02382 fftw_plan backwardPlan= 02383 fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(), 02384 (fftw_complex *)result.begin(), 02385 FFTW_BACKWARD, FFTW_ESTIMATE ); 02386 typedef typename 02387 NumericTraits<typename DestImage::Accessor::value_type>::isScalar 02388 isScalarResult; 02389 02390 // convolve with filters in freq. domain 02391 for (unsigned int i= 0; i < filters.size(); i++) 02392 { 02393 combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]), 02394 destImage(result), std::multiplies<FFTWComplex<> >()); 02395 02396 // FFT back into spatial domain (inplace in destImage) 02397 fftw_execute(backwardPlan); 02398 02399 // normalization (after FFTs), maybe stripping imaginary part 02400 applyFourierFilterImplNormalization(result, 02401 results[i].upperLeft(), results[i].accessor(), 02402 isScalarResult()); 02403 } 02404 fftw_destroy_plan(backwardPlan); 02405 } 02406 02407 /********************************************************/ 02408 /* */ 02409 /* fourierTransformReal */ 02410 /* */ 02411 /********************************************************/ 02412 02413 /** \brief Real Fourier transforms for even and odd boundary conditions 02414 (aka. cosine and sine transforms). 02415 02416 02417 If the image is real and has even symmetry, its Fourier transform 02418 is also real and has even symmetry. The Fourier transform of a real image with odd 02419 symmetry is imaginary and has odd symmetry. In either case, only about a quarter 02420 of the pixels need to be stored because the rest can be calculated from the symmetry 02421 properties. This is especially useful, if the original image is implicitly assumed 02422 to have reflective or anti-reflective boundary conditions. Then the "negative" 02423 pixel locations are defined as 02424 02425 \code 02426 even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1) 02427 odd (anti-reflective boundary conditions): f[-1] = 0 02428 f[-x] = -f[x-2] (x = 2,...,N-1) 02429 \endcode 02430 02431 end similar at the other boundary (see the FFTW documentation for details). 02432 This has the advantage that more efficient Fourier transforms that use only 02433 real numbers can be implemented. These are also known as cosine and sine transforms 02434 respectively. 02435 02436 If you use the odd transform it is important to note that in the Fourier domain, 02437 the DC component is always zero and is therefore dropped from the data structure. 02438 This means that index 0 in an odd symmetric Fourier domain image refers to 02439 the <i>first</i> harmonic. This is especially important if an image is first 02440 cosine transformed (even symmetry), then in the Fourier domain multiplied 02441 with an odd symmetric filter (e.g. a first derivative) and finally transformed 02442 back to the spatial domain with a sine transform (odd symmetric). For this to work 02443 properly the image must be shifted left or up by one pixel (depending on whether 02444 the x- or y-axis is odd symmetric) before the inverse transform can be applied. 02445 (see example below). 02446 02447 The real Fourier transform functions are named <tt>fourierTransformReal??</tt> 02448 where the questions marks stand for either <tt>E</tt> or <tt>O</tt> indicating 02449 whether the x- and y-axis is to be transformed using even or odd symmetry. 02450 The same functions can be used for both the forward and inverse transforms, 02451 only the normalization changes. For signal processing, the following 02452 normalization factors are most appropriate: 02453 02454 \code 02455 forward inverse 02456 ------------------------------------------------------------ 02457 X even, Y even 1.0 4.0 * (w-1) * (h-1) 02458 X even, Y odd -1.0 -4.0 * (w-1) * (h+1) 02459 X odd, Y even -1.0 -4.0 * (w+1) * (h-1) 02460 X odd, Y odd 1.0 4.0 * (w+1) * (h+1) 02461 \endcode 02462 02463 where <tt>w</tt> and <tt>h</tt> denote the image width and height. 02464 02465 <b> Declarations:</b> 02466 02467 pass arguments explicitly: 02468 \code 02469 namespace vigra { 02470 template <class SrcTraverser, class SrcAccessor, 02471 class DestTraverser, class DestAccessor> 02472 void 02473 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02474 DestTraverser dul, DestAccessor dest, fftw_real norm); 02475 02476 fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise 02477 } 02478 \endcode 02479 02480 02481 use argument objects in conjunction with \ref ArgumentObjectFactories : 02482 \code 02483 namespace vigra { 02484 template <class SrcTraverser, class SrcAccessor, 02485 class DestTraverser, class DestAccessor> 02486 void 02487 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02488 pair<DestTraverser, DestAccessor> dest, fftw_real norm); 02489 02490 fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise 02491 } 02492 \endcode 02493 02494 <b> Usage:</b> 02495 02496 <b>\#include</b> <vigra/fftw3.hxx><br> 02497 Namespace: vigra 02498 02499 \code 02500 vigra::FImage spatial(width,height), fourier(width,height); 02501 ... // fill image with data 02502 02503 // forward cosine transform == reflective boundary conditions 02504 fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0); 02505 02506 // multiply with a first derivative of Gaussian in x-direction 02507 for(int y = 0; y < height; ++y) 02508 { 02509 for(int x = 1; x < width; ++x) 02510 { 02511 double dx = x * M_PI / (width - 1); 02512 double dy = y * M_PI / (height - 1); 02513 fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0); 02514 } 02515 fourier(width-1, y) = 0.0; 02516 } 02517 02518 // inverse transform -- odd symmetry in x-direction, even in y, 02519 // due to symmetry of the filter 02520 fourierTransformRealOE(srcImageRange(fourier), destImage(spatial), 02521 (fftw_real)-4.0 * (width+1) * (height-1)); 02522 \endcode 02523 */ 02524 doxygen_overloaded_function(template <...> void fourierTransformReal) 02525 02526 template <class SrcTraverser, class SrcAccessor, 02527 class DestTraverser, class DestAccessor> 02528 inline void 02529 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02530 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 02531 { 02532 fourierTransformRealEE(src.first, src.second, src.third, 02533 dest.first, dest.second, norm); 02534 } 02535 02536 template <class SrcTraverser, class SrcAccessor, 02537 class DestTraverser, class DestAccessor> 02538 inline void 02539 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02540 DestTraverser dul, DestAccessor dest, fftw_real norm) 02541 { 02542 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02543 norm, FFTW_REDFT00, FFTW_REDFT00); 02544 } 02545 02546 template <class DestTraverser, class DestAccessor> 02547 inline void 02548 fourierTransformRealEE( 02549 FFTWRealImage::const_traverser sul, 02550 FFTWRealImage::const_traverser slr, 02551 FFTWRealImage::Accessor src, 02552 DestTraverser dul, DestAccessor dest, fftw_real norm) 02553 { 02554 int w = slr.x - sul.x; 02555 02556 // test for right memory layout (fftw expects a width*height fftw_real array) 02557 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 02558 fourierTransformRealImpl(sul, slr, dul, dest, 02559 norm, FFTW_REDFT00, FFTW_REDFT00); 02560 else 02561 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02562 norm, FFTW_REDFT00, FFTW_REDFT00); 02563 } 02564 02565 /********************************************************************/ 02566 02567 template <class SrcTraverser, class SrcAccessor, 02568 class DestTraverser, class DestAccessor> 02569 inline void 02570 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02571 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 02572 { 02573 fourierTransformRealOE(src.first, src.second, src.third, 02574 dest.first, dest.second, norm); 02575 } 02576 02577 template <class SrcTraverser, class SrcAccessor, 02578 class DestTraverser, class DestAccessor> 02579 inline void 02580 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02581 DestTraverser dul, DestAccessor dest, fftw_real norm) 02582 { 02583 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02584 norm, FFTW_RODFT00, FFTW_REDFT00); 02585 } 02586 02587 template <class DestTraverser, class DestAccessor> 02588 inline void 02589 fourierTransformRealOE( 02590 FFTWRealImage::const_traverser sul, 02591 FFTWRealImage::const_traverser slr, 02592 FFTWRealImage::Accessor src, 02593 DestTraverser dul, DestAccessor dest, fftw_real norm) 02594 { 02595 int w = slr.x - sul.x; 02596 02597 // test for right memory layout (fftw expects a width*height fftw_real array) 02598 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 02599 fourierTransformRealImpl(sul, slr, dul, dest, 02600 norm, FFTW_RODFT00, FFTW_REDFT00); 02601 else 02602 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02603 norm, FFTW_RODFT00, FFTW_REDFT00); 02604 } 02605 02606 /********************************************************************/ 02607 02608 template <class SrcTraverser, class SrcAccessor, 02609 class DestTraverser, class DestAccessor> 02610 inline void 02611 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02612 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 02613 { 02614 fourierTransformRealEO(src.first, src.second, src.third, 02615 dest.first, dest.second, norm); 02616 } 02617 02618 template <class SrcTraverser, class SrcAccessor, 02619 class DestTraverser, class DestAccessor> 02620 inline void 02621 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02622 DestTraverser dul, DestAccessor dest, fftw_real norm) 02623 { 02624 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02625 norm, FFTW_REDFT00, FFTW_RODFT00); 02626 } 02627 02628 template <class DestTraverser, class DestAccessor> 02629 inline void 02630 fourierTransformRealEO( 02631 FFTWRealImage::const_traverser sul, 02632 FFTWRealImage::const_traverser slr, 02633 FFTWRealImage::Accessor src, 02634 DestTraverser dul, DestAccessor dest, fftw_real norm) 02635 { 02636 int w = slr.x - sul.x; 02637 02638 // test for right memory layout (fftw expects a width*height fftw_real array) 02639 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 02640 fourierTransformRealImpl(sul, slr, dul, dest, 02641 norm, FFTW_REDFT00, FFTW_RODFT00); 02642 else 02643 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02644 norm, FFTW_REDFT00, FFTW_RODFT00); 02645 } 02646 02647 /********************************************************************/ 02648 02649 template <class SrcTraverser, class SrcAccessor, 02650 class DestTraverser, class DestAccessor> 02651 inline void 02652 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 02653 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 02654 { 02655 fourierTransformRealOO(src.first, src.second, src.third, 02656 dest.first, dest.second, norm); 02657 } 02658 02659 template <class SrcTraverser, class SrcAccessor, 02660 class DestTraverser, class DestAccessor> 02661 inline void 02662 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02663 DestTraverser dul, DestAccessor dest, fftw_real norm) 02664 { 02665 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02666 norm, FFTW_RODFT00, FFTW_RODFT00); 02667 } 02668 02669 template <class DestTraverser, class DestAccessor> 02670 inline void 02671 fourierTransformRealOO( 02672 FFTWRealImage::const_traverser sul, 02673 FFTWRealImage::const_traverser slr, 02674 FFTWRealImage::Accessor src, 02675 DestTraverser dul, DestAccessor dest, fftw_real norm) 02676 { 02677 int w = slr.x - sul.x; 02678 02679 // test for right memory layout (fftw expects a width*height fftw_real array) 02680 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1)))) 02681 fourierTransformRealImpl(sul, slr, dul, dest, 02682 norm, FFTW_RODFT00, FFTW_RODFT00); 02683 else 02684 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest, 02685 norm, FFTW_RODFT00, FFTW_RODFT00); 02686 } 02687 02688 /*******************************************************************/ 02689 02690 template <class SrcTraverser, class SrcAccessor, 02691 class DestTraverser, class DestAccessor> 02692 void 02693 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 02694 DestTraverser dul, DestAccessor dest, 02695 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy) 02696 { 02697 FFTWRealImage workImage(slr - sul); 02698 copyImage(srcIterRange(sul, slr, src), destImage(workImage)); 02699 FFTWRealImage const & cworkImage = workImage; 02700 fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), 02701 dul, dest, norm, kindx, kindy); 02702 } 02703 02704 02705 template <class DestTraverser, class DestAccessor> 02706 void 02707 fourierTransformRealImpl( 02708 FFTWRealImage::const_traverser sul, 02709 FFTWRealImage::const_traverser slr, 02710 DestTraverser dul, DestAccessor dest, 02711 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy) 02712 { 02713 int w = slr.x - sul.x; 02714 int h = slr.y - sul.y; 02715 BasicImage<fftw_real> res(w, h); 02716 02717 fftw_plan plan = fftw_plan_r2r_2d(h, w, 02718 (fftw_real *)&(*sul), (fftw_real *)res.begin(), 02719 kindy, kindx, FFTW_ESTIMATE); 02720 fftw_execute(plan); 02721 fftw_destroy_plan(plan); 02722 02723 if(norm != 1.0) 02724 transformImage(srcImageRange(res), destIter(dul, dest), 02725 std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm)); 02726 else 02727 copyImage(srcImageRange(res), destIter(dul, dest)); 02728 } 02729 02730 02731 //@} 02732 02733 } // namespace vigra 02734 02735 #endif // VIGRA_FFTW3_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|