[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/fftw.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.3.3, Aug 18 2005 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 #ifndef VIGRA_FFTW_HXX 00024 #define VIGRA_FFTW_HXX 00025 00026 #include <cmath> 00027 #include <functional> 00028 #include "vigra/stdimage.hxx" 00029 #include "vigra/copyimage.hxx" 00030 #include "vigra/transformimage.hxx" 00031 #include "vigra/combineimages.hxx" 00032 #include "vigra/numerictraits.hxx" 00033 #include "vigra/imagecontainer.hxx" 00034 #include <fftw.h> 00035 00036 namespace vigra { 00037 00038 /********************************************************/ 00039 /* */ 00040 /* FFTWComplex */ 00041 /* */ 00042 /********************************************************/ 00043 00044 /* documentation: see fftw3.hxx 00045 */ 00046 class FFTWComplex 00047 : public fftw_complex 00048 { 00049 public: 00050 /** The complex' component type, as defined in '<TT>fftw.h</TT>' 00051 */ 00052 typedef fftw_real value_type; 00053 00054 /** reference type (result of operator[]) 00055 */ 00056 typedef fftw_real & reference; 00057 00058 /** const reference type (result of operator[] const) 00059 */ 00060 typedef fftw_real const & const_reference; 00061 00062 /** iterator type (result of begin() ) 00063 */ 00064 typedef fftw_real * iterator; 00065 00066 /** const iterator type (result of begin() const) 00067 */ 00068 typedef fftw_real const * const_iterator; 00069 00070 /** The norm type (result of magnitde()) 00071 */ 00072 typedef fftw_real NormType; 00073 00074 /** The squared norm type (result of squaredMagnitde()) 00075 */ 00076 typedef fftw_real SquaredNormType; 00077 00078 /** Construct from real and imaginary part. 00079 Default: 0. 00080 */ 00081 FFTWComplex(value_type const & ire = 0.0, value_type const & iim = 0.0) 00082 { 00083 re = ire; 00084 im = iim; 00085 } 00086 00087 /** Copy constructor. 00088 */ 00089 FFTWComplex(FFTWComplex const & o) 00090 : fftw_complex(o) 00091 {} 00092 00093 /** Construct from plain <TT>fftw_complex</TT>. 00094 */ 00095 FFTWComplex(fftw_complex const & o) 00096 : fftw_complex(o) 00097 {} 00098 00099 /** Construct from TinyVector. 00100 */ 00101 template <class T> 00102 FFTWComplex(TinyVector<T, 2> const & o) 00103 { 00104 re = o[0]; 00105 im = o[1]; 00106 } 00107 00108 /** Assignment. 00109 */ 00110 FFTWComplex& operator=(FFTWComplex const & o) 00111 { 00112 re = o.re; 00113 im = o.im; 00114 return *this; 00115 } 00116 00117 /** Assignment. 00118 */ 00119 FFTWComplex& operator=(fftw_complex const & o) 00120 { 00121 re = o.re; 00122 im = o.im; 00123 return *this; 00124 } 00125 00126 /** Assignment. 00127 */ 00128 FFTWComplex& operator=(fftw_real const & o) 00129 { 00130 re = o; 00131 im = 0.0; 00132 return *this; 00133 } 00134 00135 /** Assignment. 00136 */ 00137 template <class T> 00138 FFTWComplex& operator=(TinyVector<T, 2> const & o) 00139 { 00140 re = o[0]; 00141 im = o[1]; 00142 return *this; 00143 } 00144 00145 /** Unary negation. 00146 */ 00147 FFTWComplex operator-() const 00148 { return FFTWComplex(-re, -im); } 00149 00150 /** Squared magnitude x*conj(x) 00151 */ 00152 SquaredNormType squaredMagnitude() const 00153 { return c_re(*this)*c_re(*this)+c_im(*this)*c_im(*this); } 00154 00155 /** Magnitude (length of radius vector). 00156 */ 00157 NormType magnitude() const 00158 { return VIGRA_CSTD::sqrt(squaredMagnitude()); } 00159 00160 /** Phase angle. 00161 */ 00162 value_type phase() const 00163 { return VIGRA_CSTD::atan2(c_im(*this),c_re(*this)); } 00164 00165 /** Access components as if number were a vector. 00166 */ 00167 reference operator[](int i) 00168 { return (&re)[i]; } 00169 00170 /** Read components as if number were a vector. 00171 */ 00172 const_reference operator[](int i) const 00173 { return (&re)[i]; } 00174 00175 /** Length of complex number (always 2). 00176 */ 00177 int size() const 00178 { return 2; } 00179 00180 iterator begin() 00181 { return &re; } 00182 00183 iterator end() 00184 { return &re + 2; } 00185 00186 const_iterator begin() const 00187 { return &re; } 00188 00189 const_iterator end() const 00190 { return &re + 2; } 00191 }; 00192 00193 /********************************************************/ 00194 /* */ 00195 /* FFTWComplex Traits */ 00196 /* */ 00197 /********************************************************/ 00198 00199 /* documentation: see fftw3.hxx 00200 */ 00201 template<> 00202 struct NumericTraits<fftw_complex> 00203 { 00204 typedef fftw_complex Type; 00205 typedef fftw_complex Promote; 00206 typedef fftw_complex RealPromote; 00207 typedef fftw_complex ComplexPromote; 00208 typedef fftw_real ValueType; 00209 00210 typedef VigraFalseType isIntegral; 00211 typedef VigraFalseType isScalar; 00212 typedef VigraFalseType isOrdered; 00213 typedef VigraTrueType isComplex; 00214 00215 static FFTWComplex zero() { return FFTWComplex(0.0, 0.0); } 00216 static FFTWComplex one() { return FFTWComplex(1.0, 0.0); } 00217 static FFTWComplex nonZero() { return one(); } 00218 00219 static const Promote & toPromote(const Type & v) { return v; } 00220 static const RealPromote & toRealPromote(const Type & v) { return v; } 00221 static const Type & fromPromote(const Promote & v) { return v; } 00222 static const Type & fromRealPromote(const RealPromote & v) { return v; } 00223 }; 00224 00225 template<> 00226 struct NumericTraits<FFTWComplex> 00227 : public NumericTraits<fftw_complex> 00228 { 00229 typedef FFTWComplex Type; 00230 typedef FFTWComplex Promote; 00231 typedef FFTWComplex RealPromote; 00232 typedef FFTWComplex ComplexPromote; 00233 typedef fftw_real ValueType; 00234 00235 typedef VigraFalseType isIntegral; 00236 typedef VigraFalseType isScalar; 00237 typedef VigraFalseType isOrdered; 00238 typedef VigraTrueType isComplex; 00239 }; 00240 00241 template<> 00242 struct NormTraits<fftw_complex> 00243 { 00244 typedef fftw_complex Type; 00245 typedef fftw_real SquaredNormType; 00246 typedef fftw_real NormType; 00247 }; 00248 00249 template<> 00250 struct NormTraits<FFTWComplex> 00251 { 00252 typedef FFTWComplex Type; 00253 typedef fftw_real SquaredNormType; 00254 typedef fftw_real NormType; 00255 }; 00256 00257 00258 template <> 00259 struct PromoteTraits<fftw_complex, fftw_complex> 00260 { 00261 typedef fftw_complex Promote; 00262 }; 00263 00264 template <> 00265 struct PromoteTraits<fftw_complex, double> 00266 { 00267 typedef fftw_complex Promote; 00268 }; 00269 00270 template <> 00271 struct PromoteTraits<double, fftw_complex> 00272 { 00273 typedef fftw_complex Promote; 00274 }; 00275 00276 template <> 00277 struct PromoteTraits<FFTWComplex, FFTWComplex> 00278 { 00279 typedef FFTWComplex Promote; 00280 }; 00281 00282 template <> 00283 struct PromoteTraits<FFTWComplex, double> 00284 { 00285 typedef FFTWComplex Promote; 00286 }; 00287 00288 template <> 00289 struct PromoteTraits<double, FFTWComplex> 00290 { 00291 typedef FFTWComplex Promote; 00292 }; 00293 00294 00295 /********************************************************/ 00296 /* */ 00297 /* FFTWComplex Operations */ 00298 /* */ 00299 /********************************************************/ 00300 00301 /* documentation: see fftw3.hxx 00302 */ 00303 inline bool operator ==(FFTWComplex const &a, const FFTWComplex &b) { 00304 return c_re(a) == c_re(b) && c_im(a) == c_im(b); 00305 } 00306 00307 inline bool operator !=(FFTWComplex const &a, const FFTWComplex &b) { 00308 return c_re(a) != c_re(b) || c_im(a) != c_im(b); 00309 } 00310 00311 inline FFTWComplex &operator +=(FFTWComplex &a, const FFTWComplex &b) { 00312 c_re(a) += c_re(b); 00313 c_im(a) += c_im(b); 00314 return a; 00315 } 00316 00317 inline FFTWComplex &operator -=(FFTWComplex &a, const FFTWComplex &b) { 00318 c_re(a) -= c_re(b); 00319 c_im(a) -= c_im(b); 00320 return a; 00321 } 00322 00323 inline FFTWComplex &operator *=(FFTWComplex &a, const FFTWComplex &b) { 00324 FFTWComplex::value_type t = c_re(a)*c_re(b)-c_im(a)*c_im(b); 00325 c_im(a) = c_re(a)*c_im(b)+c_im(a)*c_re(b); 00326 c_re(a) = t; 00327 return a; 00328 } 00329 00330 inline FFTWComplex &operator /=(FFTWComplex &a, const FFTWComplex &b) { 00331 FFTWComplex::value_type sm = b.squaredMagnitude(); 00332 FFTWComplex::value_type t = (c_re(a)*c_re(b)+c_im(a)*c_im(b))/sm; 00333 c_im(a) = (c_re(b)*c_im(a)-c_re(a)*c_im(b))/sm; 00334 c_re(a) = t; 00335 return a; 00336 } 00337 00338 inline FFTWComplex &operator *=(FFTWComplex &a, const double &b) { 00339 c_re(a) *= b; 00340 c_im(a) *= b; 00341 return a; 00342 } 00343 00344 inline FFTWComplex &operator /=(FFTWComplex &a, const double &b) { 00345 c_re(a) /= b; 00346 c_im(a) /= b; 00347 return a; 00348 } 00349 00350 inline FFTWComplex operator +(FFTWComplex a, const FFTWComplex &b) { 00351 a += b; 00352 return a; 00353 } 00354 00355 inline FFTWComplex operator -(FFTWComplex a, const FFTWComplex &b) { 00356 a -= b; 00357 return a; 00358 } 00359 00360 inline FFTWComplex operator *(FFTWComplex a, const FFTWComplex &b) { 00361 a *= b; 00362 return a; 00363 } 00364 00365 inline FFTWComplex operator *(FFTWComplex a, const double &b) { 00366 a *= b; 00367 return a; 00368 } 00369 00370 inline FFTWComplex operator *(const double &a, FFTWComplex b) { 00371 b *= a; 00372 return b; 00373 } 00374 00375 inline FFTWComplex operator /(FFTWComplex a, const FFTWComplex &b) { 00376 a /= b; 00377 return a; 00378 } 00379 00380 inline FFTWComplex operator /(FFTWComplex a, const double &b) { 00381 a /= b; 00382 return a; 00383 } 00384 00385 using VIGRA_CSTD::abs; 00386 00387 inline FFTWComplex::value_type abs(const FFTWComplex &a) 00388 { 00389 return a.magnitude(); 00390 } 00391 00392 inline FFTWComplex conj(const FFTWComplex &a) 00393 { 00394 return FFTWComplex(a.re, -a.im); 00395 } 00396 00397 inline FFTWComplex::NormType norm(const FFTWComplex &a) 00398 { 00399 return a.magnitude(); 00400 } 00401 00402 inline FFTWComplex::SquaredNormType squaredNorm(const FFTWComplex &a) 00403 { 00404 return a.squaredMagnitude(); 00405 } 00406 00407 /********************************************************/ 00408 /* */ 00409 /* FFTWRealImage */ 00410 /* */ 00411 /********************************************************/ 00412 00413 /* documentation: see fftw3.hxx 00414 */ 00415 typedef BasicImage<fftw_real> FFTWRealImage; 00416 00417 /********************************************************/ 00418 /* */ 00419 /* FFTWComplexImage */ 00420 /* */ 00421 /********************************************************/ 00422 00423 template<> 00424 struct IteratorTraits< 00425 BasicImageIterator<FFTWComplex, FFTWComplex **> > 00426 { 00427 typedef BasicImageIterator<FFTWComplex, FFTWComplex **> Iterator; 00428 typedef Iterator iterator; 00429 typedef iterator::iterator_category iterator_category; 00430 typedef iterator::value_type value_type; 00431 typedef iterator::reference reference; 00432 typedef iterator::index_reference index_reference; 00433 typedef iterator::pointer pointer; 00434 typedef iterator::difference_type difference_type; 00435 typedef iterator::row_iterator row_iterator; 00436 typedef iterator::column_iterator column_iterator; 00437 typedef VectorAccessor<FFTWComplex> default_accessor; 00438 typedef VectorAccessor<FFTWComplex> DefaultAccessor; 00439 typedef VigraTrueType hasConstantStrides; 00440 }; 00441 00442 template<> 00443 struct IteratorTraits< 00444 ConstBasicImageIterator<FFTWComplex, FFTWComplex **> > 00445 { 00446 typedef ConstBasicImageIterator<FFTWComplex, FFTWComplex **> Iterator; 00447 typedef Iterator iterator; 00448 typedef iterator::iterator_category iterator_category; 00449 typedef iterator::value_type value_type; 00450 typedef iterator::reference reference; 00451 typedef iterator::index_reference index_reference; 00452 typedef iterator::pointer pointer; 00453 typedef iterator::difference_type difference_type; 00454 typedef iterator::row_iterator row_iterator; 00455 typedef iterator::column_iterator column_iterator; 00456 typedef VectorAccessor<FFTWComplex> default_accessor; 00457 typedef VectorAccessor<FFTWComplex> DefaultAccessor; 00458 typedef VigraTrueType hasConstantStrides; 00459 }; 00460 00461 /* documentation: see fftw3.hxx 00462 */ 00463 typedef BasicImage<FFTWComplex> FFTWComplexImage; 00464 00465 /********************************************************/ 00466 /* */ 00467 /* FFTWComplex-Accessors */ 00468 /* */ 00469 /********************************************************/ 00470 00471 /* documentation: see fftw3.hxx 00472 */ 00473 class FFTWRealAccessor 00474 { 00475 public: 00476 00477 /// The accessor's value type. 00478 typedef fftw_real value_type; 00479 00480 /// Read real part at iterator position. 00481 template <class ITERATOR> 00482 value_type operator()(ITERATOR const & i) const { 00483 return c_re(*i); 00484 } 00485 00486 /// Read real part at offset from iterator position. 00487 template <class ITERATOR, class DIFFERENCE> 00488 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00489 return c_re(i[d]); 00490 } 00491 00492 /// Write real part at iterator position. 00493 template <class ITERATOR> 00494 void set(value_type const & v, ITERATOR const & i) const { 00495 c_re(*i)= v; 00496 } 00497 00498 /// Write real part at offset from iterator position. 00499 template <class ITERATOR, class DIFFERENCE> 00500 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 00501 c_re(i[d])= v; 00502 } 00503 }; 00504 00505 /* documentation: see fftw3.hxx 00506 */ 00507 class FFTWImaginaryAccessor 00508 { 00509 public: 00510 /// The accessor's value type. 00511 typedef fftw_real value_type; 00512 00513 /// Read imaginary part at iterator position. 00514 template <class ITERATOR> 00515 value_type operator()(ITERATOR const & i) const { 00516 return c_im(*i); 00517 } 00518 00519 /// Read imaginary part at offset from iterator position. 00520 template <class ITERATOR, class DIFFERENCE> 00521 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00522 return c_im(i[d]); 00523 } 00524 00525 /// Write imaginary part at iterator position. 00526 template <class ITERATOR> 00527 void set(value_type const & v, ITERATOR const & i) const { 00528 c_im(*i)= v; 00529 } 00530 00531 /// Write imaginary part at offset from iterator position. 00532 template <class ITERATOR, class DIFFERENCE> 00533 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 00534 c_im(i[d])= v; 00535 } 00536 }; 00537 00538 /* documentation: see fftw3.hxx 00539 */ 00540 class FFTWWriteRealAccessor: public FFTWRealAccessor 00541 { 00542 public: 00543 /// The accessor's value type. 00544 typedef fftw_real value_type; 00545 00546 /** Write real number at iterator position. Set imaginary part 00547 to 0. 00548 */ 00549 template <class ITERATOR> 00550 void set(value_type const & v, ITERATOR const & i) const { 00551 c_re(*i)= v; 00552 c_im(*i)= 0; 00553 } 00554 00555 /** Write real number at offset from iterator position. Set imaginary part 00556 to 0. 00557 */ 00558 template <class ITERATOR, class DIFFERENCE> 00559 void set(value_type const & v, ITERATOR const & i, DIFFERENCE d) const { 00560 c_re(i[d])= v; 00561 c_im(i[d])= 0; 00562 } 00563 }; 00564 00565 /* documentation: see fftw3.hxx 00566 */ 00567 class FFTWMagnitudeAccessor 00568 { 00569 public: 00570 /// The accessor's value type. 00571 typedef fftw_real value_type; 00572 00573 /// Read magnitude at iterator position. 00574 template <class ITERATOR> 00575 value_type operator()(ITERATOR const & i) const { 00576 return (*i).magnitude(); 00577 } 00578 00579 /// Read magnitude at offset from iterator position. 00580 template <class ITERATOR, class DIFFERENCE> 00581 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00582 return (i[d]).magnitude(); 00583 } 00584 }; 00585 00586 /* documentation: see fftw3.hxx 00587 */ 00588 class FFTWPhaseAccessor 00589 { 00590 public: 00591 /// The accessor's value type. 00592 typedef fftw_real value_type; 00593 00594 /// Read phase at iterator position. 00595 template <class ITERATOR> 00596 value_type operator()(ITERATOR const & i) const { 00597 return (*i).phase(); 00598 } 00599 00600 /// Read phase at offset from iterator position. 00601 template <class ITERATOR, class DIFFERENCE> 00602 value_type operator()(ITERATOR const & i, DIFFERENCE d) const { 00603 return (i[d]).phase(); 00604 } 00605 }; 00606 00607 /********************************************************/ 00608 /* */ 00609 /* Fourier Transform */ 00610 /* */ 00611 /********************************************************/ 00612 00613 /** \page FourierTransformFFTW2 Fast Fourier Transform 00614 00615 This documentation describes the deprecated VIGRA interface to 00616 FFTW 2. Use the \link FourierTransform interface to the newer 00617 version FFTW 3\endlink instead. 00618 00619 VIGRA uses the <a href="http://www.fftw.org/">FFTW Fast Fourier 00620 Transform</a> package to perform Fourier transformations. VIGRA 00621 provides a wrapper for FFTW's complex number type (FFTWComplex), 00622 but FFTW's functions are used verbatim. If the image is stored as 00623 a FFTWComplexImage, a FFT is performed like this: 00624 00625 \code 00626 vigra::FFTWComplexImage spatial(width,height), fourier(width,height); 00627 ... // fill image with data 00628 00629 // create a plan for optimal performance 00630 fftwnd_plan forwardPlan= 00631 fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE ); 00632 00633 // calculate FFT 00634 fftwnd_one(forwardPlan, spatial.begin(), fourier.begin()); 00635 \endcode 00636 00637 Note that in the creation of a plan, the height must be given 00638 first. Note also that <TT>spatial.begin()</TT> may only be passed 00639 to <TT>fftwnd_one</TT> if the transform shall be applied to the 00640 entire image. When you want to retrict operation to an ROI, you 00641 create a copy of the ROI in an image of appropriate size. 00642 00643 More information on using FFTW can be found <a href="http://www.fftw.org/doc/">here</a>. 00644 00645 FFTW produces fourier images that have the DC component (the 00646 origin of the Fourier space) in the upper left corner. Often, one 00647 wants the origin in the center of the image, so that frequencies 00648 always increase towards the border of the image. This can be 00649 achieved by calling \ref moveDCToCenter(). The inverse 00650 transformation is done by \ref moveDCToUpperLeft(). 00651 00652 <b>\#include</b> "<a href="fftw_8hxx-source.html">vigra/fftw.hxx</a>"<br> 00653 Namespace: vigra 00654 */ 00655 00656 /********************************************************/ 00657 /* */ 00658 /* moveDCToCenter */ 00659 /* */ 00660 /********************************************************/ 00661 00662 /* documentation: see fftw3.hxx 00663 */ 00664 template <class SrcImageIterator, class SrcAccessor, 00665 class DestImageIterator, class DestAccessor> 00666 void moveDCToCenter(SrcImageIterator src_upperleft, 00667 SrcImageIterator src_lowerright, SrcAccessor sa, 00668 DestImageIterator dest_upperleft, DestAccessor da) 00669 { 00670 int w= src_lowerright.x - src_upperleft.x; 00671 int h= src_lowerright.y - src_upperleft.y; 00672 int w1 = w/2; 00673 int h1 = h/2; 00674 int w2 = (w+1)/2; 00675 int h2 = (h+1)/2; 00676 00677 // 2. Quadrant zum 4. 00678 copyImage(srcIterRange(src_upperleft, 00679 src_upperleft + Diff2D(w2, h2), sa), 00680 destIter (dest_upperleft + Diff2D(w1, h1), da)); 00681 00682 // 4. Quadrant zum 2. 00683 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2), 00684 src_lowerright, sa), 00685 destIter (dest_upperleft, da)); 00686 00687 // 1. Quadrant zum 3. 00688 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0), 00689 src_upperleft + Diff2D(w, h2), sa), 00690 destIter (dest_upperleft + Diff2D(0, h1), da)); 00691 00692 // 3. Quadrant zum 1. 00693 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2), 00694 src_upperleft + Diff2D(w2, h), sa), 00695 destIter (dest_upperleft + Diff2D(w1, 0), da)); 00696 } 00697 00698 template <class SrcImageIterator, class SrcAccessor, 00699 class DestImageIterator, class DestAccessor> 00700 inline void moveDCToCenter( 00701 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00702 pair<DestImageIterator, DestAccessor> dest) 00703 { 00704 moveDCToCenter(src.first, src.second, src.third, 00705 dest.first, dest.second); 00706 } 00707 00708 /********************************************************/ 00709 /* */ 00710 /* moveDCToUpperLeft */ 00711 /* */ 00712 /********************************************************/ 00713 00714 /* documentation: see fftw3.hxx 00715 */ 00716 template <class SrcImageIterator, class SrcAccessor, 00717 class DestImageIterator, class DestAccessor> 00718 void moveDCToUpperLeft(SrcImageIterator src_upperleft, 00719 SrcImageIterator src_lowerright, SrcAccessor sa, 00720 DestImageIterator dest_upperleft, DestAccessor da) 00721 { 00722 int w= src_lowerright.x - src_upperleft.x; 00723 int h= src_lowerright.y - src_upperleft.y; 00724 int w2 = w/2; 00725 int h2 = h/2; 00726 int w1 = (w+1)/2; 00727 int h1 = (h+1)/2; 00728 00729 // 2. Quadrant zum 4. 00730 copyImage(srcIterRange(src_upperleft, 00731 src_upperleft + Diff2D(w2, h2), sa), 00732 destIter (dest_upperleft + Diff2D(w1, h1), da)); 00733 00734 // 4. Quadrant zum 2. 00735 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2), 00736 src_lowerright, sa), 00737 destIter (dest_upperleft, da)); 00738 00739 // 1. Quadrant zum 3. 00740 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0), 00741 src_upperleft + Diff2D(w, h2), sa), 00742 destIter (dest_upperleft + Diff2D(0, h1), da)); 00743 00744 // 3. Quadrant zum 1. 00745 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2), 00746 src_upperleft + Diff2D(w2, h), sa), 00747 destIter (dest_upperleft + Diff2D(w1, 0), da)); 00748 } 00749 00750 template <class SrcImageIterator, class SrcAccessor, 00751 class DestImageIterator, class DestAccessor> 00752 inline void moveDCToUpperLeft( 00753 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00754 pair<DestImageIterator, DestAccessor> dest) 00755 { 00756 moveDCToUpperLeft(src.first, src.second, src.third, 00757 dest.first, dest.second); 00758 } 00759 00760 /********************************************************/ 00761 /* */ 00762 /* applyFourierFilter */ 00763 /* */ 00764 /********************************************************/ 00765 00766 /* documentation: see fftw3.hxx 00767 */ 00768 00769 // applyFourierFilter versions without fftwnd_plans: 00770 template <class SrcImageIterator, class SrcAccessor, 00771 class FilterImageIterator, class FilterAccessor, 00772 class DestImageIterator, class DestAccessor> 00773 inline 00774 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00775 pair<FilterImageIterator, FilterAccessor> filter, 00776 pair<DestImageIterator, DestAccessor> dest) 00777 { 00778 applyFourierFilter(src.first, src.second, src.third, 00779 filter.first, filter.second, 00780 dest.first, dest.second); 00781 } 00782 00783 template <class SrcImageIterator, class SrcAccessor, 00784 class FilterImageIterator, class FilterAccessor, 00785 class DestImageIterator, class DestAccessor> 00786 void applyFourierFilter(SrcImageIterator srcUpperLeft, 00787 SrcImageIterator srcLowerRight, SrcAccessor sa, 00788 FilterImageIterator filterUpperLeft, FilterAccessor fa, 00789 DestImageIterator destUpperLeft, DestAccessor da) 00790 { 00791 // copy real input images into a complex one... 00792 int w= srcLowerRight.x - srcUpperLeft.x; 00793 int h= srcLowerRight.y - srcUpperLeft.y; 00794 00795 FFTWComplexImage workImage(w, h); 00796 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 00797 destImage(workImage, FFTWWriteRealAccessor())); 00798 00799 // ...and call the impl 00800 FFTWComplexImage const & cworkImage = workImage; 00801 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 00802 filterUpperLeft, fa, 00803 destUpperLeft, da); 00804 } 00805 00806 typedef FFTWComplexImage::const_traverser FFTWConstTraverser; 00807 00808 template <class FilterImageIterator, class FilterAccessor, 00809 class DestImageIterator, class DestAccessor> 00810 inline 00811 void applyFourierFilter( 00812 FFTWComplexImage::const_traverser srcUpperLeft, 00813 FFTWComplexImage::const_traverser srcLowerRight, 00814 FFTWComplexImage::ConstAccessor sa, 00815 FilterImageIterator filterUpperLeft, FilterAccessor fa, 00816 DestImageIterator destUpperLeft, DestAccessor da) 00817 { 00818 int w= srcLowerRight.x - srcUpperLeft.x; 00819 int h= srcLowerRight.y - srcUpperLeft.y; 00820 00821 // test for right memory layout (fftw expects a 2*width*height floats array) 00822 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 00823 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa, 00824 filterUpperLeft, fa, 00825 destUpperLeft, da); 00826 else 00827 { 00828 FFTWComplexImage workImage(w, h); 00829 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 00830 destImage(workImage)); 00831 00832 FFTWComplexImage const & cworkImage = workImage; 00833 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 00834 filterUpperLeft, fa, 00835 destUpperLeft, da); 00836 } 00837 } 00838 00839 template <class FilterImageIterator, class FilterAccessor, 00840 class DestImageIterator, class DestAccessor> 00841 void applyFourierFilterImpl( 00842 FFTWComplexImage::const_traverser srcUpperLeft, 00843 FFTWComplexImage::const_traverser srcLowerRight, 00844 FFTWComplexImage::ConstAccessor sa, 00845 FilterImageIterator filterUpperLeft, FilterAccessor fa, 00846 DestImageIterator destUpperLeft, DestAccessor da) 00847 { 00848 // create plans and call variant with plan parameters 00849 int w= srcLowerRight.x - srcUpperLeft.x; 00850 int h= srcLowerRight.y - srcUpperLeft.y; 00851 00852 fftwnd_plan forwardPlan= 00853 fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE ); 00854 fftwnd_plan backwardPlan= 00855 fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE); 00856 00857 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa, 00858 filterUpperLeft, fa, 00859 destUpperLeft, da, 00860 forwardPlan, backwardPlan); 00861 00862 fftwnd_destroy_plan(forwardPlan); 00863 fftwnd_destroy_plan(backwardPlan); 00864 } 00865 00866 // applyFourierFilter versions with fftwnd_plans: 00867 template <class SrcImageIterator, class SrcAccessor, 00868 class FilterImageIterator, class FilterAccessor, 00869 class DestImageIterator, class DestAccessor> 00870 inline 00871 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 00872 pair<FilterImageIterator, FilterAccessor> filter, 00873 pair<DestImageIterator, DestAccessor> dest, 00874 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 00875 { 00876 applyFourierFilter(src.first, src.second, src.third, 00877 filter.first, filter.second, 00878 dest.first, dest.second, 00879 forwardPlan, backwardPlan); 00880 } 00881 00882 template <class SrcImageIterator, class SrcAccessor, 00883 class FilterImageIterator, class FilterAccessor, 00884 class DestImageIterator, class DestAccessor> 00885 void applyFourierFilter(SrcImageIterator srcUpperLeft, 00886 SrcImageIterator srcLowerRight, SrcAccessor sa, 00887 FilterImageIterator filterUpperLeft, FilterAccessor fa, 00888 DestImageIterator destUpperLeft, DestAccessor da, 00889 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 00890 { 00891 int w= srcLowerRight.x - srcUpperLeft.x; 00892 int h= srcLowerRight.y - srcUpperLeft.y; 00893 00894 FFTWComplexImage workImage(w, h); 00895 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 00896 destImage(workImage, FFTWWriteRealAccessor())); 00897 00898 FFTWComplexImage const & cworkImage = workImage; 00899 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 00900 filterUpperLeft, fa, 00901 destUpperLeft, da, 00902 forwardPlan, backwardPlan); 00903 } 00904 00905 template <class FilterImageIterator, class FilterAccessor, 00906 class DestImageIterator, class DestAccessor> 00907 inline 00908 void applyFourierFilter( 00909 FFTWComplexImage::const_traverser srcUpperLeft, 00910 FFTWComplexImage::const_traverser srcLowerRight, 00911 FFTWComplexImage::ConstAccessor sa, 00912 FilterImageIterator filterUpperLeft, FilterAccessor fa, 00913 DestImageIterator destUpperLeft, DestAccessor da, 00914 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 00915 { 00916 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa, 00917 filterUpperLeft, fa, 00918 destUpperLeft, da, 00919 forwardPlan, backwardPlan); 00920 } 00921 00922 template <class FilterImageIterator, class FilterAccessor, 00923 class DestImageIterator, class DestAccessor> 00924 void applyFourierFilterImpl( 00925 FFTWComplexImage::const_traverser srcUpperLeft, 00926 FFTWComplexImage::const_traverser srcLowerRight, 00927 FFTWComplexImage::ConstAccessor sa, 00928 FilterImageIterator filterUpperLeft, FilterAccessor fa, 00929 DestImageIterator destUpperLeft, DestAccessor da, 00930 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 00931 { 00932 FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft); 00933 00934 // FFT from srcImage to complexResultImg 00935 fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)), 00936 complexResultImg.begin()); 00937 00938 // convolve in freq. domain (in complexResultImg) 00939 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa), 00940 destImage(complexResultImg), std::multiplies<FFTWComplex>()); 00941 00942 // FFT back into spatial domain (inplace in complexResultImg) 00943 fftwnd_one(backwardPlan, complexResultImg.begin(), 0); 00944 00945 typedef typename 00946 NumericTraits<typename DestAccessor::value_type>::isScalar 00947 isScalarResult; 00948 00949 // normalization (after FFTs), maybe stripping imaginary part 00950 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da, 00951 isScalarResult()); 00952 } 00953 00954 template <class DestImageIterator, class DestAccessor> 00955 void applyFourierFilterImplNormalization(FFTWComplexImage const &srcImage, 00956 DestImageIterator destUpperLeft, 00957 DestAccessor da, 00958 VigraFalseType) 00959 { 00960 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 00961 00962 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 00963 { 00964 DestImageIterator dIt= destUpperLeft; 00965 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 00966 { 00967 da.setComponent(srcImage(x, y).re*normFactor, dIt, 0); 00968 da.setComponent(srcImage(x, y).im*normFactor, dIt, 1); 00969 } 00970 } 00971 } 00972 00973 inline 00974 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 00975 FFTWComplexImage::traverser destUpperLeft, 00976 FFTWComplexImage::Accessor da, 00977 VigraFalseType) 00978 { 00979 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da), 00980 linearIntensityTransform<FFTWComplex>(1.0/(srcImage.width() * srcImage.height()))); 00981 } 00982 00983 template <class DestImageIterator, class DestAccessor> 00984 void applyFourierFilterImplNormalization(FFTWComplexImage const & srcImage, 00985 DestImageIterator destUpperLeft, 00986 DestAccessor da, 00987 VigraTrueType) 00988 { 00989 double normFactor= 1.0/(srcImage.width() * srcImage.height()); 00990 00991 for(int y=0; y<srcImage.height(); y++, destUpperLeft.y++) 00992 { 00993 DestImageIterator dIt= destUpperLeft; 00994 for(int x= 0; x< srcImage.width(); x++, dIt.x++) 00995 da.set(srcImage(x, y).re*normFactor, dIt); 00996 } 00997 } 00998 00999 /**********************************************************/ 01000 /* */ 01001 /* applyFourierFilterFamily */ 01002 /* */ 01003 /**********************************************************/ 01004 01005 /* documentation: see fftw3.hxx 01006 */ 01007 01008 // applyFourierFilterFamily versions without fftwnd_plans: 01009 template <class SrcImageIterator, class SrcAccessor, 01010 class FilterType, class DestImage> 01011 inline 01012 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01013 const ImageArray<FilterType> &filters, 01014 ImageArray<DestImage> &results) 01015 { 01016 applyFourierFilterFamily(src.first, src.second, src.third, 01017 filters, results); 01018 } 01019 01020 template <class SrcImageIterator, class SrcAccessor, 01021 class FilterType, class DestImage> 01022 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 01023 SrcImageIterator srcLowerRight, SrcAccessor sa, 01024 const ImageArray<FilterType> &filters, 01025 ImageArray<DestImage> &results) 01026 { 01027 int w= srcLowerRight.x - srcUpperLeft.x; 01028 int h= srcLowerRight.y - srcUpperLeft.y; 01029 01030 FFTWComplexImage workImage(w, h); 01031 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01032 destImage(workImage, FFTWWriteRealAccessor())); 01033 01034 FFTWComplexImage const & cworkImage = workImage; 01035 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01036 filters, results); 01037 } 01038 01039 template <class FilterType, class DestImage> 01040 inline 01041 void applyFourierFilterFamily( 01042 FFTWComplexImage::const_traverser srcUpperLeft, 01043 FFTWComplexImage::const_traverser srcLowerRight, 01044 FFTWComplexImage::ConstAccessor sa, 01045 const ImageArray<FilterType> &filters, 01046 ImageArray<DestImage> &results) 01047 { 01048 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa, 01049 filters, results); 01050 } 01051 01052 template <class FilterType, class DestImage> 01053 void applyFourierFilterFamilyImpl( 01054 FFTWComplexImage::const_traverser srcUpperLeft, 01055 FFTWComplexImage::const_traverser srcLowerRight, 01056 FFTWComplexImage::ConstAccessor sa, 01057 const ImageArray<FilterType> &filters, 01058 ImageArray<DestImage> &results) 01059 { 01060 int w= srcLowerRight.x - srcUpperLeft.x; 01061 int h= srcLowerRight.y - srcUpperLeft.y; 01062 01063 fftwnd_plan forwardPlan= 01064 fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_ESTIMATE ); 01065 fftwnd_plan backwardPlan= 01066 fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE); 01067 01068 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa, 01069 filters, results, 01070 forwardPlan, backwardPlan); 01071 01072 fftwnd_destroy_plan(forwardPlan); 01073 fftwnd_destroy_plan(backwardPlan); 01074 } 01075 01076 // applyFourierFilterFamily versions with fftwnd_plans: 01077 template <class SrcImageIterator, class SrcAccessor, 01078 class FilterType, class DestImage> 01079 inline 01080 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, 01081 const ImageArray<FilterType> &filters, 01082 ImageArray<DestImage> &results, 01083 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 01084 { 01085 applyFourierFilterFamily(src.first, src.second, src.third, 01086 filters, results, 01087 forwardPlan, backwardPlan); 01088 } 01089 01090 template <class SrcImageIterator, class SrcAccessor, 01091 class FilterType, class DestImage> 01092 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft, 01093 SrcImageIterator srcLowerRight, SrcAccessor sa, 01094 const ImageArray<FilterType> &filters, 01095 ImageArray<DestImage> &results, 01096 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 01097 { 01098 int w= srcLowerRight.x - srcUpperLeft.x; 01099 int h= srcLowerRight.y - srcUpperLeft.y; 01100 01101 FFTWComplexImage workImage(w, h); 01102 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01103 destImage(workImage, FFTWWriteRealAccessor())); 01104 01105 FFTWComplexImage const & cworkImage = workImage; 01106 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01107 filters, results, 01108 forwardPlan, backwardPlan); 01109 } 01110 01111 template <class FilterType, class DestImage> 01112 inline 01113 void applyFourierFilterFamily( 01114 FFTWComplexImage::const_traverser srcUpperLeft, 01115 FFTWComplexImage::const_traverser srcLowerRight, 01116 FFTWComplexImage::ConstAccessor sa, 01117 const ImageArray<FilterType> &filters, 01118 ImageArray<DestImage> &results, 01119 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 01120 { 01121 int w= srcLowerRight.x - srcUpperLeft.x; 01122 int h= srcLowerRight.y - srcUpperLeft.y; 01123 01124 // test for right memory layout (fftw expects a 2*width*height floats array) 01125 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1)))) 01126 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa, 01127 filters, results, 01128 forwardPlan, backwardPlan); 01129 else 01130 { 01131 FFTWComplexImage workImage(w, h); 01132 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa), 01133 destImage(workImage)); 01134 01135 FFTWComplexImage const & cworkImage = workImage; 01136 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(), 01137 filters, results, 01138 forwardPlan, backwardPlan); 01139 } 01140 } 01141 01142 template <class FilterType, class DestImage> 01143 void applyFourierFilterFamilyImpl( 01144 FFTWComplexImage::const_traverser srcUpperLeft, 01145 FFTWComplexImage::const_traverser srcLowerRight, 01146 FFTWComplexImage::ConstAccessor sa, 01147 const ImageArray<FilterType> &filters, 01148 ImageArray<DestImage> &results, 01149 const fftwnd_plan &forwardPlan, const fftwnd_plan &backwardPlan) 01150 { 01151 // make sure the filter images have the right dimensions 01152 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(), 01153 "applyFourierFilterFamily called with src image size != filters.imageSize()!"); 01154 01155 // make sure the result image array has the right dimensions 01156 results.resize(filters.size()); 01157 results.resizeImages(filters.imageSize()); 01158 01159 // FFT from srcImage to freqImage 01160 int w= srcLowerRight.x - srcUpperLeft.x; 01161 int h= srcLowerRight.y - srcUpperLeft.y; 01162 01163 FFTWComplexImage freqImage(w, h); 01164 FFTWComplexImage result(w, h); 01165 01166 fftwnd_one(forwardPlan, const_cast<FFTWComplex *>(&(*srcUpperLeft)), freqImage.begin()); 01167 01168 typedef typename 01169 NumericTraits<typename DestImage::Accessor::value_type>::isScalar 01170 isScalarResult; 01171 01172 // convolve with filters in freq. domain 01173 for (unsigned int i= 0; i < filters.size(); i++) 01174 { 01175 combineTwoImages(srcImageRange(freqImage), srcImage(filters[i]), 01176 destImage(result), std::multiplies<FFTWComplex>()); 01177 01178 // FFT back into spatial domain (inplace in destImage) 01179 fftwnd_one(backwardPlan, result.begin(), 0); 01180 01181 // normalization (after FFTs), maybe stripping imaginary part 01182 applyFourierFilterImplNormalization(result, 01183 results[i].upperLeft(), results[i].accessor(), 01184 isScalarResult()); 01185 } 01186 } 01187 01188 } // namespace vigra 01189 01190 #endif // VIGRA_FFTW_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|