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