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