36 #ifndef VIGRA_FFTW3_HXX 37 #define VIGRA_FFTW3_HXX 42 #include "stdimage.hxx" 43 #include "copyimage.hxx" 44 #include "transformimage.hxx" 45 #include "combineimages.hxx" 46 #include "numerictraits.hxx" 47 #include "imagecontainer.hxx" 52 typedef double fftw_real;
58 struct FFTWReal<fftw_complex>
64 struct FFTWReal<fftwf_complex>
70 struct FFTWReal<fftwl_complex>
72 typedef long double type;
76 struct FFTWReal2Complex;
79 struct FFTWReal2Complex<double>
81 typedef fftw_complex type;
82 typedef fftw_plan plan_type;
86 struct FFTWReal2Complex<float>
88 typedef fftwf_complex type;
89 typedef fftwf_plan plan_type;
93 struct FFTWReal2Complex<long double>
95 typedef fftwl_complex type;
96 typedef fftwl_plan plan_type;
130 template <
class Real =
double>
169 FFTWComplex(value_type
const & re = 0.0, value_type
const & im = 0.0)
179 data_[0] = o.data_[0];
180 data_[1] = o.data_[1];
188 data_[0] = (Real)o.real();
189 data_[1] = (Real)o.imag();
196 data_[0] = (Real)o[0];
197 data_[1] = (Real)o[1];
204 data_[0] = (Real)o[0];
205 data_[1] = (Real)o[1];
212 data_[0] = (Real)o[0];
213 data_[1] = (Real)o[1];
221 data_[0] = (Real)o.real();
222 data_[1] = (Real)o.imag();
230 data_[0] = (Real)o[0];
231 data_[1] = (Real)o[1];
238 data_[0] = o.data_[0];
239 data_[1] = o.data_[1];
248 data_[0] = (Real)o.real();
249 data_[1] = (Real)o.imag();
257 data_[0] = (Real)o[0];
258 data_[1] = (Real)o[1];
266 data_[0] = (Real)o[0];
267 data_[1] = (Real)o[1];
275 data_[0] = (Real)o[0];
276 data_[1] = (Real)o[1];
312 data_[0] = (Real)o[0];
313 data_[1] = (Real)o[1];
322 data_[0] = (Real)o.real();
323 data_[1] = (Real)o.imag();
330 const_reference re()
const 336 const_reference
real()
const 342 const_reference im()
const 348 const_reference
imag()
const 359 {
return data_[0]*data_[0]+data_[1]*data_[1]; }
390 {
return data_ + 2; }
392 const_iterator begin()
const 395 const_iterator end()
const 396 {
return data_ + 2; }
508 struct NumericTraits<fftw_complex>
510 typedef fftw_complex Type;
511 typedef fftw_complex Promote;
512 typedef fftw_complex RealPromote;
513 typedef fftw_complex ComplexPromote;
514 typedef fftw_real ValueType;
516 typedef VigraFalseType isIntegral;
517 typedef VigraFalseType isScalar;
518 typedef NumericTraits<fftw_real>::isSigned isSigned;
519 typedef VigraFalseType isOrdered;
520 typedef VigraTrueType isComplex;
526 static const Promote & toPromote(
const Type & v) {
return v; }
527 static const RealPromote & toRealPromote(
const Type & v) {
return v; }
528 static const Type & fromPromote(
const Promote & v) {
return v; }
529 static const Type & fromRealPromote(
const RealPromote & v) {
return v; }
533 struct NumericTraits<FFTWComplex<Real> >
539 typedef typename Type::value_type ValueType;
541 typedef VigraFalseType isIntegral;
542 typedef VigraFalseType isScalar;
543 typedef typename NumericTraits<ValueType>::isSigned isSigned;
544 typedef VigraFalseType isOrdered;
545 typedef VigraTrueType isComplex;
547 static Type zero() {
return Type(0.0, 0.0); }
548 static Type one() {
return Type(1.0, 0.0); }
549 static Type nonZero() {
return one(); }
551 static const Promote & toPromote(
const Type & v) {
return v; }
552 static const RealPromote & toRealPromote(
const Type & v) {
return v; }
553 static const Type & fromPromote(
const Promote & v) {
return v; }
554 static const Type & fromRealPromote(
const RealPromote & v) {
return v; }
558 struct NormTraits<fftw_complex>
560 typedef fftw_complex Type;
566 struct NormTraits<FFTWComplex<Real> >
570 typedef typename Type::NormType
NormType;
574 struct PromoteTraits<fftw_complex, fftw_complex>
576 typedef fftw_complex Promote;
580 struct PromoteTraits<fftw_complex, double>
582 typedef fftw_complex Promote;
586 struct PromoteTraits<double, fftw_complex>
588 typedef fftw_complex Promote;
591 template <
class Real>
592 struct PromoteTraits<FFTWComplex<Real>, FFTWComplex<Real> >
597 template <
class Real>
598 struct PromoteTraits<FFTWComplex<Real>, double>
603 template <
class Real>
604 struct PromoteTraits<double, FFTWComplex<Real> >
610 struct CanSkipInitialization<std::complex<T> >
612 typedef typename CanSkipInitialization<T>::type type;
613 static const bool value = type::asBool;
617 struct CanSkipInitialization<FFTWComplex<Real> >
619 typedef typename CanSkipInitialization<Real>::type type;
620 static const bool value = type::asBool;
623 namespace multi_math {
626 struct MultiMathOperand;
628 template <
class Real>
629 struct MultiMathOperand<FFTWComplex<Real> >
631 typedef MultiMathOperand<FFTWComplex<Real> > AllowOverload;
634 static const int ndim = 0;
640 template <
class SHAPE>
641 bool checkShape(SHAPE
const &)
const 646 template <
class SHAPE>
652 void inc(
unsigned int )
const 655 void reset(
unsigned int )
const 672 typedef std::size_t size_type;
673 typedef std::ptrdiff_t difference_type;
675 typedef const Ty *const_pointer;
676 typedef Ty& reference;
677 typedef const Ty& const_reference;
678 typedef Ty value_type;
680 pointer address(reference val)
const 683 const_pointer address(const_reference val)
const 686 template<
class Other>
689 typedef FFTWAllocator<Other> other;
692 FFTWAllocator()
throw()
695 template<
class Other>
696 FFTWAllocator(
const FFTWAllocator<Other>& right)
throw()
699 template<
class Other>
700 FFTWAllocator& operator=(
const FFTWAllocator<Other>& right)
705 pointer allocate(size_type count,
void * = 0)
707 return (pointer)fftw_malloc(count *
sizeof(Ty));
710 void deallocate(pointer ptr, size_type count)
715 void construct(pointer ptr,
const Ty& val)
721 void destroy(pointer ptr)
726 size_type max_size()
const throw()
728 return NumericTraits<std::ptrdiff_t>::max() /
sizeof(Ty);
741 typedef size_t size_type;
742 typedef std::ptrdiff_t difference_type;
743 typedef value_type *pointer;
744 typedef const value_type *const_pointer;
745 typedef value_type& reference;
746 typedef const value_type& const_reference;
748 pointer address(reference val)
const 751 const_pointer address(const_reference val)
const 754 template<
class Other>
757 typedef allocator<Other> other;
763 template<
class Other>
764 allocator(
const allocator<Other>& right)
throw()
767 template<
class Other>
768 allocator& operator=(
const allocator<Other>& right)
773 pointer allocate(size_type count,
void * = 0)
775 return (pointer)fftw_malloc(count *
sizeof(value_type));
778 void deallocate(pointer ptr, size_type )
783 void construct(pointer ptr,
const value_type& val)
785 new(ptr) value_type(val);
789 void destroy(pointer ptr)
794 size_type max_size()
const throw()
796 return vigra::NumericTraits<std::ptrdiff_t>::max() /
sizeof(value_type);
826 return a.re() == b.re() && a.im() == b.im();
831 return a.re() == b && a.im() == 0.0;
836 return a == b.re() && 0.0 == b.im();
842 return a.re() != b.re() || a.im() != b.im();
848 return a.re() != b || a.im() != 0.0;
854 return a != b.re() || 0.0 != b.im();
877 a.im() = a.re()*b.im()+a.im()*b.re();
887 a.im() = (b.re()*a.im()-a.re()*b.im())/sm;
1049 #define VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(fct) \ 1050 template <class R> \ 1051 inline FFTWComplex<R> fct(const FFTWComplex<R> &a) \ 1053 return std::fct(reinterpret_cast<std::complex<R> const &>(a)); \ 1056 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
cos)
1057 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(cosh)
1058 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
exp)
1059 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
log)
1060 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
log10)
1061 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
sin)
1062 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(sinh)
1063 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
sqrt)
1064 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(
tan)
1065 VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION(tanh)
1067 #undef VIGRA_DEFINE_FFTW_COMPLEX_FUNCTION 1072 return std::pow(
reinterpret_cast<std::complex<R>
const &
>(a), e);
1078 return std::pow(
reinterpret_cast<std::complex<R>
const &
>(a), e);
1084 return std::pow(
reinterpret_cast<std::complex<R>
const &
>(a),
1085 reinterpret_cast<std::complex<R>
const &
>(e));
1091 return std::pow(a,
reinterpret_cast<std::complex<R>
const &
>(e));
1100 template <
class Real>
1101 ostream & operator<<(ostream & s, vigra::FFTWComplex<Real>
const & v)
1103 s << std::complex<Real>(v.re(), v.im());
1146 typedef Iterator iterator;
1149 typedef typename iterator::iterator_category iterator_category;
1150 typedef typename iterator::value_type value_type;
1151 typedef typename iterator::reference reference;
1152 typedef typename iterator::index_reference index_reference;
1153 typedef typename iterator::pointer pointer;
1154 typedef typename iterator::difference_type difference_type;
1155 typedef typename iterator::row_iterator row_iterator;
1156 typedef typename iterator::column_iterator column_iterator;
1159 typedef VigraTrueType hasConstantStrides;
1164 ConstBasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> >
1168 typedef Iterator iterator;
1171 typedef typename iterator::iterator_category iterator_category;
1172 typedef typename iterator::value_type value_type;
1173 typedef typename iterator::reference reference;
1174 typedef typename iterator::index_reference index_reference;
1175 typedef typename iterator::pointer pointer;
1176 typedef typename iterator::difference_type difference_type;
1177 typedef typename iterator::row_iterator row_iterator;
1178 typedef typename iterator::column_iterator column_iterator;
1181 typedef VigraTrueType hasConstantStrides;
1216 template <
class Real =
double>
1225 template <
class ITERATOR>
1231 template <
class ITERATOR,
class DIFFERENCE>
1237 template <
class ITERATOR>
1238 void set(value_type
const & v, ITERATOR
const & i)
const {
1243 template <
class ITERATOR,
class DIFFERENCE>
1244 void set(value_type
const & v, ITERATOR
const & i, DIFFERENCE d)
const {
1249 template <
class R,
class ITERATOR>
1255 template <
class R,
class ITERATOR,
class DIFFERENCE>
1267 template <
class Real =
double>
1275 template <
class ITERATOR>
1281 template <
class ITERATOR,
class DIFFERENCE>
1287 template <
class ITERATOR>
1288 void set(value_type
const & v, ITERATOR
const & i)
const {
1293 template <
class ITERATOR,
class DIFFERENCE>
1294 void set(value_type
const & v, ITERATOR
const & i, DIFFERENCE d)
const {
1299 template <
class R,
class ITERATOR>
1305 template <
class R,
class ITERATOR,
class DIFFERENCE>
1318 template <
class Real =
double>
1329 template <
class ITERATOR>
1330 void set(value_type
const & v, ITERATOR
const & i)
const {
1338 template <
class ITERATOR,
class DIFFERENCE>
1339 void set(value_type
const & v, ITERATOR
const & i, DIFFERENCE d)
const {
1350 template <
class Real =
double>
1358 template <
class ITERATOR>
1360 return (*i).squaredMagnitude();
1364 template <
class ITERATOR,
class DIFFERENCE>
1366 return (i[d]).squaredMagnitude();
1376 template <
class Real =
double>
1384 template <
class ITERATOR>
1386 return (*i).magnitude();
1390 template <
class ITERATOR,
class DIFFERENCE>
1392 return (i[d]).magnitude();
1402 template <
class Real =
double>
1410 template <
class ITERATOR>
1412 return (*i).phase();
1416 template <
class ITERATOR,
class DIFFERENCE>
1418 return (i[d]).phase();
1594 template <
class SrcImageIterator,
class SrcAccessor,
1595 class DestImageIterator,
class DestAccessor>
1596 void moveDCToCenter(SrcImageIterator src_upperleft,
1597 SrcImageIterator src_lowerright, SrcAccessor sa,
1598 DestImageIterator dest_upperleft, DestAccessor da)
1600 int w = int(src_lowerright.x - src_upperleft.x);
1601 int h = int(src_lowerright.y - src_upperleft.y);
1609 src_upperleft +
Diff2D(w2, h2), sa),
1610 destIter (dest_upperleft +
Diff2D(w1, h1), da));
1614 src_lowerright, sa),
1615 destIter (dest_upperleft, da));
1619 src_upperleft +
Diff2D(w, h2), sa),
1620 destIter (dest_upperleft +
Diff2D(0, h1), da));
1624 src_upperleft +
Diff2D(w2, h), sa),
1625 destIter (dest_upperleft +
Diff2D(w1, 0), da));
1628 template <
class SrcImageIterator,
class SrcAccessor,
1629 class DestImageIterator,
class DestAccessor>
1630 inline void moveDCToCenter(
1631 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1632 pair<DestImageIterator, DestAccessor> dest)
1634 moveDCToCenter(src.first, src.second, src.third,
1635 dest.first, dest.second);
1685 template <
class SrcImageIterator,
class SrcAccessor,
1686 class DestImageIterator,
class DestAccessor>
1687 void moveDCToUpperLeft(SrcImageIterator src_upperleft,
1688 SrcImageIterator src_lowerright, SrcAccessor sa,
1689 DestImageIterator dest_upperleft, DestAccessor da)
1691 int w = int(src_lowerright.x - src_upperleft.x);
1692 int h = int(src_lowerright.y - src_upperleft.y);
1700 src_upperleft +
Diff2D(w2, h2), sa),
1701 destIter (dest_upperleft +
Diff2D(w1, h1), da));
1705 src_lowerright, sa),
1706 destIter (dest_upperleft, da));
1710 src_upperleft +
Diff2D(w, h2), sa),
1711 destIter (dest_upperleft +
Diff2D(0, h1), da));
1715 src_upperleft +
Diff2D(w2, h), sa),
1716 destIter (dest_upperleft +
Diff2D(w1, 0), da));
1719 template <
class SrcImageIterator,
class SrcAccessor,
1720 class DestImageIterator,
class DestAccessor>
1721 inline void moveDCToUpperLeft(
1722 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1723 pair<DestImageIterator, DestAccessor> dest)
1725 moveDCToUpperLeft(src.first, src.second, src.third,
1726 dest.first, dest.second);
1737 int w = int(slr.x - sul.x);
1738 int h = int(slr.y - sul.y);
1740 FFTWComplexImage sworkImage, dworkImage;
1742 fftw_complex * srcPtr = (fftw_complex *)(&*sul);
1743 fftw_complex * destPtr = (fftw_complex *)(&*dul);
1746 if (h > 1 && &(*(sul +
Diff2D(w, 0))) != &(*(sul +
Diff2D(0, 1))))
1749 copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
1750 srcPtr = (fftw_complex *)(&(*sworkImage.
upperLeft()));
1752 if (h > 1 && &(*(dul +
Diff2D(w, 0))) != &(*(dul +
Diff2D(0, 1))))
1755 destPtr = (fftw_complex *)(&(*dworkImage.
upperLeft()));
1758 fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
1760 fftw_destroy_plan(plan);
1762 if (h > 1 && &(*(dul +
Diff2D(w, 0))) != &(*(dul +
Diff2D(0, 1))))
1764 copyImage(srcImageRange(dworkImage), destIter(dul, dest));
1923 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
1926 template <
class SrcImageIterator,
class SrcAccessor>
1927 void fourierTransform(SrcImageIterator srcUpperLeft,
1928 SrcImageIterator srcLowerRight, SrcAccessor sa,
1932 int w= srcLowerRight.x - srcUpperLeft.x;
1933 int h= srcLowerRight.y - srcUpperLeft.y;
1935 FFTWComplexImage workImage(w, h);
1936 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1940 FFTWComplexImage
const & cworkImage = workImage;
1945 template <
class SrcImageIterator,
class SrcAccessor>
1947 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1948 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
1950 fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
1964 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
1967 template <
class DestImageIterator,
class DestAccessor>
1970 DestImageIterator dul, DestAccessor dest)
1972 int w = slr.x - sul.x;
1973 int h = slr.y - sul.y;
1975 FFTWComplexImage workImage(w, h);
1976 fourierTransformInverse(sul, slr, src, workImage.
upperLeft(), workImage.
accessor());
1977 copyImage(srcImageRange(workImage), destIter(dul, dest));
1981 template <
class DestImageIterator,
class DestAccessor>
1985 pair<DestImageIterator, DestAccessor> dest)
1987 fourierTransformInverse(src.first, src.second, src.third, dest.first, dest.second);
2088 template <
class SrcImageIterator,
class SrcAccessor,
2089 class FilterImageIterator,
class FilterAccessor,
2090 class DestImageIterator,
class DestAccessor>
2091 void applyFourierFilter(SrcImageIterator srcUpperLeft,
2092 SrcImageIterator srcLowerRight, SrcAccessor sa,
2093 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2094 DestImageIterator destUpperLeft, DestAccessor da)
2097 int w = int(srcLowerRight.x - srcUpperLeft.x);
2098 int h = int(srcLowerRight.y - srcUpperLeft.y);
2100 FFTWComplexImage workImage(w, h);
2101 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2105 FFTWComplexImage
const & cworkImage = workImage;
2107 filterUpperLeft, fa,
2111 template <
class FilterImageIterator,
class FilterAccessor,
2112 class DestImageIterator,
class DestAccessor>
2114 void applyFourierFilter(
2118 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2119 DestImageIterator destUpperLeft, DestAccessor da)
2121 int w = srcLowerRight.x - srcUpperLeft.x;
2122 int h = srcLowerRight.y - srcUpperLeft.y;
2125 if (&(*(srcUpperLeft +
Diff2D(w, 0))) == &(*(srcUpperLeft +
Diff2D(0, 1))))
2126 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
2127 filterUpperLeft, fa,
2131 FFTWComplexImage workImage(w, h);
2132 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2133 destImage(workImage));
2135 FFTWComplexImage
const & cworkImage = workImage;
2137 filterUpperLeft, fa,
2142 template <
class SrcImageIterator,
class SrcAccessor,
2143 class FilterImageIterator,
class FilterAccessor,
2144 class DestImageIterator,
class DestAccessor>
2146 void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2147 pair<FilterImageIterator, FilterAccessor> filter,
2148 pair<DestImageIterator, DestAccessor> dest)
2150 applyFourierFilter(src.first, src.second, src.third,
2151 filter.first, filter.second,
2152 dest.first, dest.second);
2155 template <
class FilterImageIterator,
class FilterAccessor,
2156 class DestImageIterator,
class DestAccessor>
2157 void applyFourierFilterImpl(
2161 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2162 DestImageIterator destUpperLeft, DestAccessor da)
2164 int w = int(srcLowerRight.x - srcUpperLeft.x);
2165 int h = int(srcLowerRight.y - srcUpperLeft.y);
2167 FFTWComplexImage complexResultImg(srcLowerRight - srcUpperLeft);
2170 fftw_plan forwardPlan=
2171 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2172 (fftw_complex *)complexResultImg.
begin(),
2173 FFTW_FORWARD, FFTW_ESTIMATE );
2174 fftw_execute(forwardPlan);
2175 fftw_destroy_plan(forwardPlan);
2178 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
2179 destImage(complexResultImg), std::multiplies<
FFTWComplex<> >());
2182 fftw_plan backwardPlan=
2183 fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.
begin(),
2184 (fftw_complex *)complexResultImg.
begin(),
2185 FFTW_BACKWARD, FFTW_ESTIMATE);
2186 fftw_execute(backwardPlan);
2187 fftw_destroy_plan(backwardPlan);
2190 NumericTraits<typename DestAccessor::value_type>::isScalar
2194 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
2198 template <
class DestImageIterator,
class DestAccessor>
2199 void applyFourierFilterImplNormalization(FFTWComplexImage
const &srcImage,
2200 DestImageIterator destUpperLeft,
2204 double normFactor= 1.0/(srcImage.
width() * srcImage.
height());
2206 for(
int y=0; y<srcImage.
height(); y++, destUpperLeft.y++)
2208 DestImageIterator dIt= destUpperLeft;
2209 for(
int x= 0; x< srcImage.
width(); x++, dIt.x++)
2211 da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
2212 da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
2218 void applyFourierFilterImplNormalization(FFTWComplexImage
const & srcImage,
2223 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
2227 template <
class DestImageIterator,
class DestAccessor>
2228 void applyFourierFilterImplNormalization(FFTWComplexImage
const & srcImage,
2229 DestImageIterator destUpperLeft,
2233 double normFactor= 1.0/(srcImage.
width() * srcImage.
height());
2235 for(
int y=0; y<srcImage.
height(); y++, destUpperLeft.y++)
2237 DestImageIterator dIt= destUpperLeft;
2238 for(
int x= 0; x< srcImage.
width(); x++, dIt.x++)
2239 da.set(srcImage(x, y).re()*normFactor, dIt);
2315 template <
class SrcImageIterator,
class SrcAccessor,
2316 class FilterType,
class DestImage>
2318 void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
2322 applyFourierFilterFamily(src.first, src.second, src.third,
2326 template <
class SrcImageIterator,
class SrcAccessor,
2327 class FilterType,
class DestImage>
2328 void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
2329 SrcImageIterator srcLowerRight, SrcAccessor sa,
2333 int w = int(srcLowerRight.x - srcUpperLeft.x);
2334 int h = int(srcLowerRight.y - srcUpperLeft.y);
2336 FFTWComplexImage workImage(w, h);
2337 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2340 FFTWComplexImage
const & cworkImage = workImage;
2345 template <
class FilterType,
class DestImage>
2347 void applyFourierFilterFamily(
2354 int w= srcLowerRight.x - srcUpperLeft.x;
2357 if (&(*(srcUpperLeft +
Diff2D(w, 0))) == &(*(srcUpperLeft +
Diff2D(0, 1))))
2358 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
2362 int h = srcLowerRight.y - srcUpperLeft.y;
2363 FFTWComplexImage workImage(w, h);
2364 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2365 destImage(workImage));
2367 FFTWComplexImage
const & cworkImage = workImage;
2373 template <
class FilterType,
class DestImage>
2374 void applyFourierFilterFamilyImpl(
2385 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.
imageSize(),
2386 "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2393 int w = int(srcLowerRight.x - srcUpperLeft.x);
2394 int h = int(srcLowerRight.y - srcUpperLeft.y);
2396 FFTWComplexImage freqImage(w, h);
2397 FFTWComplexImage result(w, h);
2399 fftw_plan forwardPlan=
2400 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2401 (fftw_complex *)freqImage.
begin(),
2402 FFTW_FORWARD, FFTW_ESTIMATE );
2403 fftw_execute(forwardPlan);
2404 fftw_destroy_plan(forwardPlan);
2406 fftw_plan backwardPlan=
2407 fftw_plan_dft_2d(h, w, (fftw_complex *)result.
begin(),
2408 (fftw_complex *)result.
begin(),
2409 FFTW_BACKWARD, FFTW_ESTIMATE );
2411 NumericTraits<typename DestImage::Accessor::value_type>::isScalar
2415 for (
unsigned int i= 0; i < filters.
size(); i++)
2421 fftw_execute(backwardPlan);
2424 applyFourierFilterImplNormalization(result,
2425 results[i].upperLeft(), results[i].accessor(),
2428 fftw_destroy_plan(backwardPlan);
2563 template <
class SrcTraverser,
class SrcAccessor,
2564 class DestTraverser,
class DestAccessor>
2566 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2567 pair<DestTraverser, DestAccessor> dest, fftw_real
norm)
2569 fourierTransformRealEE(src.first, src.second, src.third,
2570 dest.first, dest.second, norm);
2573 template <
class SrcTraverser,
class SrcAccessor,
2574 class DestTraverser,
class DestAccessor>
2576 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2577 DestTraverser dul, DestAccessor dest, fftw_real norm)
2579 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2580 norm, FFTW_REDFT00, FFTW_REDFT00);
2583 template <
class DestTraverser,
class DestAccessor>
2585 fourierTransformRealEE(
2589 DestTraverser dul, DestAccessor dest, fftw_real norm)
2591 int w = slr.x - sul.x;
2594 if (&(*(sul +
Diff2D(w, 0))) == &(*(sul +
Diff2D(0, 1))))
2595 fourierTransformRealImpl(sul, slr, dul, dest,
2596 norm, FFTW_REDFT00, FFTW_REDFT00);
2598 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2599 norm, FFTW_REDFT00, FFTW_REDFT00);
2604 template <
class SrcTraverser,
class SrcAccessor,
2605 class DestTraverser,
class DestAccessor>
2607 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2608 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2610 fourierTransformRealOE(src.first, src.second, src.third,
2611 dest.first, dest.second, norm);
2614 template <
class SrcTraverser,
class SrcAccessor,
2615 class DestTraverser,
class DestAccessor>
2617 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2618 DestTraverser dul, DestAccessor dest, fftw_real norm)
2620 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2621 norm, FFTW_RODFT00, FFTW_REDFT00);
2624 template <
class DestTraverser,
class DestAccessor>
2626 fourierTransformRealOE(
2630 DestTraverser dul, DestAccessor dest, fftw_real norm)
2632 int w = slr.x - sul.x;
2635 if (&(*(sul +
Diff2D(w, 0))) == &(*(sul +
Diff2D(0, 1))))
2636 fourierTransformRealImpl(sul, slr, dul, dest,
2637 norm, FFTW_RODFT00, FFTW_REDFT00);
2639 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2640 norm, FFTW_RODFT00, FFTW_REDFT00);
2645 template <
class SrcTraverser,
class SrcAccessor,
2646 class DestTraverser,
class DestAccessor>
2648 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2649 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2651 fourierTransformRealEO(src.first, src.second, src.third,
2652 dest.first, dest.second, norm);
2655 template <
class SrcTraverser,
class SrcAccessor,
2656 class DestTraverser,
class DestAccessor>
2658 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2659 DestTraverser dul, DestAccessor dest, fftw_real norm)
2661 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2662 norm, FFTW_REDFT00, FFTW_RODFT00);
2665 template <
class DestTraverser,
class DestAccessor>
2667 fourierTransformRealEO(
2671 DestTraverser dul, DestAccessor dest, fftw_real norm)
2673 int w = slr.x - sul.x;
2676 if (&(*(sul +
Diff2D(w, 0))) == &(*(sul +
Diff2D(0, 1))))
2677 fourierTransformRealImpl(sul, slr, dul, dest,
2678 norm, FFTW_REDFT00, FFTW_RODFT00);
2680 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2681 norm, FFTW_REDFT00, FFTW_RODFT00);
2686 template <
class SrcTraverser,
class SrcAccessor,
2687 class DestTraverser,
class DestAccessor>
2689 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2690 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2692 fourierTransformRealOO(src.first, src.second, src.third,
2693 dest.first, dest.second, norm);
2696 template <
class SrcTraverser,
class SrcAccessor,
2697 class DestTraverser,
class DestAccessor>
2699 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2700 DestTraverser dul, DestAccessor dest, fftw_real norm)
2702 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2703 norm, FFTW_RODFT00, FFTW_RODFT00);
2706 template <
class DestTraverser,
class DestAccessor>
2708 fourierTransformRealOO(
2712 DestTraverser dul, DestAccessor dest, fftw_real norm)
2714 int w = slr.x - sul.x;
2717 if (&(*(sul +
Diff2D(w, 0))) == &(*(sul +
Diff2D(0, 1))))
2718 fourierTransformRealImpl(sul, slr, dul, dest,
2719 norm, FFTW_RODFT00, FFTW_RODFT00);
2721 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2722 norm, FFTW_RODFT00, FFTW_RODFT00);
2727 template <
class SrcTraverser,
class SrcAccessor,
2728 class DestTraverser,
class DestAccessor>
2730 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2731 DestTraverser dul, DestAccessor dest,
2732 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2734 FFTWRealImage workImage(slr - sul);
2735 copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2736 FFTWRealImage
const & cworkImage = workImage;
2738 dul, dest,
norm, kindx, kindy);
2742 template <
class DestTraverser,
class DestAccessor>
2744 fourierTransformRealImpl(
2747 DestTraverser dul, DestAccessor dest,
2748 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2750 int w = slr.x - sul.x;
2751 int h = slr.y - sul.y;
2754 fftw_plan plan = fftw_plan_r2r_2d(h, w,
2755 (fftw_real *)&(*sul), (fftw_real *)res.
begin(),
2756 kindy, kindx, FFTW_ESTIMATE);
2758 fftw_destroy_plan(plan);
2762 std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2764 copyImage(srcImageRange(res), destIter(dul, dest));
2772 #endif // VIGRA_FFTW3_HXX FFTWComplex(fftwf_complex const &o)
Definition: fftw3.hxx:202
Size2D imageSize() const
Definition: imagecontainer.hxx:411
iterator begin()
Definition: basicimage.hxx:963
Definition: fftw3.hxx:1217
FFTWComplex< R > & operator*=(FFTWComplex< R > &a, double b)
multiply-assignment with scalar double
Definition: fftw3.hxx:908
FFTWComplex< R > & operator/=(FFTWComplex< R > &a, double b)
divide-assignment with scalar double
Definition: fftw3.hxx:916
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1654
BasicImage< fftw_real > FFTWRealImage
Definition: fftw3.hxx:1132
Fundamental class template for arrays of equal-sized images.
Definition: imagecontainer.hxx:72
Real value_type
Definition: fftw3.hxx:140
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read phase at offset from iterator position.
Definition: fftw3.hxx:1417
Export associated information for each image iterator.
Definition: iteratortraits.hxx:109
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition: fftw3.hxx:1030
R imag(const FFTWComplex< R > &a)
imaginary part
Definition: fftw3.hxx:1023
FFTWComplex(value_type const &re=0.0, value_type const &im=0.0)
Definition: fftw3.hxx:169
value_type SquaredNormType
Definition: fftw3.hxx:164
Accessor accessor()
Definition: basicimage.hxx:1064
R arg(const FFTWComplex< R > &a)
pahse
Definition: fftw3.hxx:1009
FFTWReal2Complex< Real >::type complex_type
Definition: fftw3.hxx:136
Definition: fftw3.hxx:1351
SquaredNormType squaredMagnitude() const
Definition: fftw3.hxx:358
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1381
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1324
Definition: fftw3.hxx:1319
value_type phase() const
Definition: fftw3.hxx:368
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
FFTWComplex & operator=(fftw_complex const &o)
Definition: fftw3.hxx:255
FFTWComplex & operator=(TinyVector< T, 2 > const &o)
Definition: fftw3.hxx:310
value_type const * const_iterator
Definition: fftw3.hxx:156
int width() const
Definition: basicimage.hxx:838
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1355
Definition: array_vector.hxx:903
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
FFTWComplex & operator=(fftwl_complex const &o)
Definition: fftw3.hxx:273
int size() const
Definition: fftw3.hxx:383
R real(const FFTWComplex< R > &a)
real part
Definition: fftw3.hxx:1016
FFTWComplex< R > & operator+=(FFTWComplex< R > &a, double b)
add-assignment with scalar double
Definition: fftw3.hxx:894
Two dimensional difference vector.
Definition: diff2d.hxx:185
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1222
FFTWComplex & operator=(std::complex< T > const &o)
Definition: fftw3.hxx:320
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1272
FFTWComplex & operator=(FFTWComplex const &o)
Definition: fftw3.hxx:236
Definition: accessor.hxx:43
value_type operator()(ITERATOR const &i) const
Read squared magnitude at iterator position.
Definition: fftw3.hxx:1359
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read squared magnitude at offset from iterator position.
Definition: fftw3.hxx:1365
Definition: fftw3.hxx:1403
int height() const
Definition: basicimage.hxx:845
Accessor for items that are STL compatible vectors.
Definition: accessor.hxx:771
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
FFTWComplex< R > & operator-=(FFTWComplex< R > &a, double b)
subtract-assignment with scalar double
Definition: fftw3.hxx:901
FFTWComplex & operator=(float o)
Definition: fftw3.hxx:291
Definition: basicimage.hxx:262
bool operator!=(double a, const FFTWComplex< R > &b)
not equal
Definition: fftw3.hxx:853
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read real part at offset from iterator position.
Definition: fftw3.hxx:1232
linalg::TemporaryMatrix< T > log10(MultiArrayView< 2, T, C > const &v)
FFTWComplex(FFTWComplex const &o)
Definition: fftw3.hxx:177
reference operator[](int i)
Definition: fftw3.hxx:373
value_type operator()(ITERATOR const &i) const
Read imaginary part at iterator position.
Definition: fftw3.hxx:1276
value_type const & const_reference
Definition: fftw3.hxx:148
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read imaginary part at offset from iterator position.
Definition: fftw3.hxx:1282
FFTWComplex(FFTWComplex< U > const &o)
Definition: fftw3.hxx:186
value_type operator()(ITERATOR const &i) const
Read phase at iterator position.
Definition: fftw3.hxx:1411
FFTWComplex & operator=(long double o)
Definition: fftw3.hxx:300
void resize(int width, int height)
Definition: basicimage.hxx:776
FFTWComplex & operator=(FFTWComplex< U > const &o)
Definition: fftw3.hxx:246
void combineTwoImages(...)
Combine two source images into destination image.
FFTWComplex(fftw_complex const &o)
Definition: fftw3.hxx:194
size_type size() const
Definition: imagecontainer.hxx:228
void copyImage(...)
Copy source image into destination image.
FFTWComplex(fftwl_complex const &o)
Definition: fftw3.hxx:210
value_type operator()(ITERATOR const &i) const
Read real part at iterator position.
Definition: fftw3.hxx:1226
FFTWComplex< R > operator*(double a, FFTWComplex< R > b)
left multiplication with scalar double
Definition: fftw3.hxx:979
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
value_type operator()(ITERATOR const &i) const
Read magnitude at iterator position.
Definition: fftw3.hxx:1385
value_type NormType
Definition: fftw3.hxx:160
Definition: fftw3.hxx:1268
traverser lowerRight()
Definition: basicimage.hxx:934
Definition: fftw3.hxx:1377
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
Fundamental class template for images.
Definition: basicimage.hxx:473
FFTWComplex operator-() const
Definition: fftw3.hxx:353
Definition: basicimage.hxx:294
FFTWComplex< R > operator-(double a, FFTWComplex< R > const &b)
left subtraction with scalar double
Definition: fftw3.hxx:959
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
void resize(size_type newSize)
Definition: imagecontainer.hxx:310
linalg::TemporaryMatrix< T > tan(MultiArrayView< 2, T, C > const &v)
value_type * iterator
Definition: fftw3.hxx:152
const_reference operator[](int i) const
Definition: fftw3.hxx:378
NormType magnitude() const
Definition: fftw3.hxx:363
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
T sign(T t)
The sign function.
Definition: mathutil.hxx:553
virtual void resizeImages(const Diff2D &newSize)
Definition: imagecontainer.hxx:418
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1407
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read magnitude at offset from iterator position.
Definition: fftw3.hxx:1391
Wrapper class for the FFTW complex types 'fftw_complex'.
Definition: fftw3.hxx:131
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
value_type & reference
Definition: fftw3.hxx:144
FFTWComplex< R > operator/(FFTWComplex< R > a, double b)
right division with scalar double
Definition: fftw3.hxx:993
FFTWComplex(std::complex< T > const &o)
Definition: fftw3.hxx:219
BasicImage< FFTWComplex<> > FFTWComplexImage
Definition: fftw3.hxx:1192
FFTWComplex(TinyVector< T, 2 > const &o)
Definition: fftw3.hxx:228
FFTWComplex & operator=(double o)
Definition: fftw3.hxx:282
traverser upperLeft()
Definition: basicimage.hxx:923
FFTWComplex< R > operator+(double a, FFTWComplex< R > b)
left addition with scalar double
Definition: fftw3.hxx:938
FFTWComplex & operator=(fftwf_complex const &o)
Definition: fftw3.hxx:264