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>
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();
359 {
return data_[0]*data_[0]+data_[1]*data_[1]; }
390 {
return data_ + 2; }
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;
522 static FFTWComplex<> zero() {
return FFTWComplex<>(0.0, 0.0); }
523 static FFTWComplex<> one() {
return FFTWComplex<>(1.0, 0.0); }
524 static FFTWComplex<> nonZero() {
return one(); }
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> >
535 typedef FFTWComplex<Real> Type;
536 typedef FFTWComplex<Real> Promote;
537 typedef FFTWComplex<Real> RealPromote;
538 typedef FFTWComplex<Real> ComplexPromote;
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;
561 typedef fftw_real SquaredNormType;
562 typedef fftw_real NormType;
566 struct NormTraits<FFTWComplex<Real> >
568 typedef FFTWComplex<Real> Type;
569 typedef typename Type::SquaredNormType SquaredNormType;
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> >
594 typedef FFTWComplex<Real> Promote;
597 template <
class Real>
598 struct PromoteTraits<FFTWComplex<Real>, double>
600 typedef FFTWComplex<Real> Promote;
603 template <
class Real>
604 struct PromoteTraits<double, FFTWComplex<Real> >
606 typedef FFTWComplex<Real> Promote;
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);
737 class allocator<vigra::FFTWComplex<Real> >
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();
830 inline bool operator ==(FFTWComplex<R>
const &a,
double b) {
831 return a.re() == b && a.im() == 0.0;
835 inline bool operator ==(
double a,
const FFTWComplex<R> &b) {
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
1070 inline FFTWComplex<R> pow(
const FFTWComplex<R> &a,
int e)
1072 return std::pow(
reinterpret_cast<std::complex<R>
const &
>(a), e);
1076 inline FFTWComplex<R> pow(
const FFTWComplex<R> &a, R
const & e)
1078 return std::pow(
reinterpret_cast<std::complex<R>
const &
>(a), e);
1082 inline FFTWComplex<R> pow(
const FFTWComplex<R> &a,
const FFTWComplex<R> & e)
1084 return std::pow(
reinterpret_cast<std::complex<R>
const &
>(a),
1085 reinterpret_cast<std::complex<R>
const &
>(e));
1089 inline FFTWComplex<R> pow(R
const & a,
const FFTWComplex<R> &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;
1163 struct IteratorTraits<
1164 ConstBasicImageIterator<FFTWComplex<R>, FFTWComplex<R> **> >
1166 typedef FFTWComplex<R> Type;
1167 typedef ConstBasicImageIterator<Type, Type **> Iterator;
1168 typedef Iterator iterator;
1169 typedef BasicImageIterator<Type, Type **> mutable_iterator;
1170 typedef ConstBasicImageIterator<Type, Type **> const_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;
1179 typedef VectorAccessor<Type> default_accessor;
1180 typedef VectorAccessor<Type> DefaultAccessor;
1181 typedef VigraTrueType hasConstantStrides;
1216 template <
class Real =
double>
1225 template <
class ITERATOR>
1231 template <
class ITERATOR,
class DIFFERENCE>
1237 template <
class ITERATOR>
1243 template <
class ITERATOR,
class DIFFERENCE>
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>
1293 template <
class ITERATOR,
class DIFFERENCE>
1299 template <
class R,
class ITERATOR>
1305 template <
class R,
class ITERATOR,
class DIFFERENCE>
1318 template <
class Real =
double>
1329 template <
class ITERATOR>
1338 template <
class ITERATOR,
class DIFFERENCE>
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>
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));
1613 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1614 src_lowerright, sa),
1615 destIter (dest_upperleft, da));
1618 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1619 src_upperleft + Diff2D(w, h2), sa),
1620 destIter (dest_upperleft + Diff2D(0, h1), da));
1623 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
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>
1631 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1632 pair<DestImageIterator, DestAccessor> dest)
1635 dest.first, dest.second);
1685 template <
class SrcImageIterator,
class SrcAccessor,
1686 class DestImageIterator,
class DestAccessor>
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));
1704 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1705 src_lowerright, sa),
1706 destIter (dest_upperleft, da));
1709 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1710 src_upperleft + Diff2D(w, h2), sa),
1711 destIter (dest_upperleft + Diff2D(0, h1), da));
1714 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
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>
1722 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1723 pair<DestImageIterator, DestAccessor> dest)
1726 dest.first, dest.second);
1737 int w = int(slr.x - sul.x);
1738 int h = int(slr.y - sul.y);
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))))
1748 sworkImage.resize(w, h);
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))))
1754 dworkImage.resize(w, h);
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>
1928 SrcImageIterator srcLowerRight, SrcAccessor sa,
1932 int w= srcLowerRight.x - srcUpperLeft.x;
1933 int h= srcLowerRight.y - srcUpperLeft.y;
1936 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1937 destImage(workImage, FFTWWriteRealAccessor<>()));
1941 fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
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;
1977 copyImage(srcImageRange(workImage), destIter(dul, dest));
1981 template <
class DestImageIterator,
class DestAccessor>
1985 pair<DestImageIterator, DestAccessor> dest)
2075 template <
class SrcImageIterator,
class SrcAccessor,
2076 class FilterImageIterator,
class FilterAccessor,
2077 class DestImageIterator,
class DestAccessor>
2079 SrcImageIterator srcLowerRight, SrcAccessor sa,
2080 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2081 DestImageIterator destUpperLeft, DestAccessor da)
2084 int w = int(srcLowerRight.x - srcUpperLeft.x);
2085 int h = int(srcLowerRight.y - srcUpperLeft.y);
2088 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2089 destImage(workImage, FFTWWriteRealAccessor<>()));
2093 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2094 filterUpperLeft, fa,
2098 template <
class FilterImageIterator,
class FilterAccessor,
2099 class DestImageIterator,
class DestAccessor>
2105 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2106 DestImageIterator destUpperLeft, DestAccessor da)
2108 int w = srcLowerRight.x - srcUpperLeft.x;
2109 int h = srcLowerRight.y - srcUpperLeft.y;
2112 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2113 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
2114 filterUpperLeft, fa,
2119 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2120 destImage(workImage));
2123 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2124 filterUpperLeft, fa,
2129 template <
class SrcImageIterator,
class SrcAccessor,
2130 class FilterImageIterator,
class FilterAccessor,
2131 class DestImageIterator,
class DestAccessor>
2134 pair<FilterImageIterator, FilterAccessor> filter,
2135 pair<DestImageIterator, DestAccessor> dest)
2138 filter.first, filter.second,
2139 dest.first, dest.second);
2142 template <
class FilterImageIterator,
class FilterAccessor,
2143 class DestImageIterator,
class DestAccessor>
2144 void applyFourierFilterImpl(
2148 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2149 DestImageIterator destUpperLeft, DestAccessor da)
2151 int w = int(srcLowerRight.x - srcUpperLeft.x);
2152 int h = int(srcLowerRight.y - srcUpperLeft.y);
2157 fftw_plan forwardPlan=
2158 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2159 (fftw_complex *)complexResultImg.begin(),
2160 FFTW_FORWARD, FFTW_ESTIMATE );
2161 fftw_execute(forwardPlan);
2162 fftw_destroy_plan(forwardPlan);
2165 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
2166 destImage(complexResultImg), std::multiplies<FFTWComplex<> >());
2169 fftw_plan backwardPlan=
2170 fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
2171 (fftw_complex *)complexResultImg.begin(),
2172 FFTW_BACKWARD, FFTW_ESTIMATE);
2173 fftw_execute(backwardPlan);
2174 fftw_destroy_plan(backwardPlan);
2177 NumericTraits<typename DestAccessor::value_type>::isScalar
2181 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
2185 template <
class DestImageIterator,
class DestAccessor>
2187 DestImageIterator destUpperLeft,
2191 double normFactor= 1.0/(srcImage.width() * srcImage.height());
2193 for(
int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2195 DestImageIterator dIt= destUpperLeft;
2196 for(
int x= 0; x< srcImage.width(); x++, dIt.x++)
2198 da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
2199 da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
2205 void applyFourierFilterImplNormalization(
FFTWComplexImage const & srcImage,
2210 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
2214 template <
class DestImageIterator,
class DestAccessor>
2215 void applyFourierFilterImplNormalization(
FFTWComplexImage const & srcImage,
2216 DestImageIterator destUpperLeft,
2220 double normFactor= 1.0/(srcImage.width() * srcImage.height());
2222 for(
int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2224 DestImageIterator dIt= destUpperLeft;
2225 for(
int x= 0; x< srcImage.width(); x++, dIt.x++)
2226 da.set(srcImage(x, y).re()*normFactor, dIt);
2291 template <
class SrcImageIterator,
class SrcAccessor,
2292 class FilterType,
class DestImage>
2295 const ImageArray<FilterType> &filters,
2296 ImageArray<DestImage> &results)
2302 template <
class SrcImageIterator,
class SrcAccessor,
2303 class FilterType,
class DestImage>
2305 SrcImageIterator srcLowerRight, SrcAccessor sa,
2306 const ImageArray<FilterType> &filters,
2307 ImageArray<DestImage> &results)
2309 int w = int(srcLowerRight.x - srcUpperLeft.x);
2310 int h = int(srcLowerRight.y - srcUpperLeft.y);
2313 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2314 destImage(workImage, FFTWWriteRealAccessor<>()));
2317 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2321 template <
class FilterType,
class DestImage>
2327 const ImageArray<FilterType> &filters,
2328 ImageArray<DestImage> &results)
2330 int w= srcLowerRight.x - srcUpperLeft.x;
2333 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2334 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
2338 int h = srcLowerRight.y - srcUpperLeft.y;
2340 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2341 destImage(workImage));
2344 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2349 template <
class FilterType,
class DestImage>
2350 void applyFourierFilterFamilyImpl(
2354 const ImageArray<FilterType> &filters,
2355 ImageArray<DestImage> &results)
2361 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
2362 "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2365 results.resize(filters.size());
2366 results.resizeImages(filters.imageSize());
2369 int w = int(srcLowerRight.x - srcUpperLeft.x);
2370 int h = int(srcLowerRight.y - srcUpperLeft.y);
2375 fftw_plan forwardPlan=
2376 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2377 (fftw_complex *)freqImage.begin(),
2378 FFTW_FORWARD, FFTW_ESTIMATE );
2379 fftw_execute(forwardPlan);
2380 fftw_destroy_plan(forwardPlan);
2382 fftw_plan backwardPlan=
2383 fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
2384 (fftw_complex *)result.begin(),
2385 FFTW_BACKWARD, FFTW_ESTIMATE );
2387 NumericTraits<typename DestImage::Accessor::value_type>::isScalar
2391 for (
unsigned int i= 0; i < filters.size(); i++)
2394 destImage(result), std::multiplies<FFTWComplex<> >());
2397 fftw_execute(backwardPlan);
2400 applyFourierFilterImplNormalization(result,
2401 results[i].upperLeft(), results[i].accessor(),
2404 fftw_destroy_plan(backwardPlan);
2526 template <
class SrcTraverser,
class SrcAccessor,
2527 class DestTraverser,
class DestAccessor>
2529 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2530 pair<DestTraverser, DestAccessor> dest, fftw_real
norm)
2532 fourierTransformRealEE(src.first, src.second, src.third,
2533 dest.first, dest.second, norm);
2536 template <
class SrcTraverser,
class SrcAccessor,
2537 class DestTraverser,
class DestAccessor>
2539 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2540 DestTraverser dul, DestAccessor dest, fftw_real
norm)
2542 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2543 norm, FFTW_REDFT00, FFTW_REDFT00);
2546 template <
class DestTraverser,
class DestAccessor>
2548 fourierTransformRealEE(
2552 DestTraverser dul, DestAccessor dest, fftw_real
norm)
2554 int w = slr.x - sul.x;
2557 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2558 fourierTransformRealImpl(sul, slr, dul, dest,
2559 norm, FFTW_REDFT00, FFTW_REDFT00);
2561 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2562 norm, FFTW_REDFT00, FFTW_REDFT00);
2567 template <
class SrcTraverser,
class SrcAccessor,
2568 class DestTraverser,
class DestAccessor>
2570 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2571 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2573 fourierTransformRealOE(src.first, src.second, src.third,
2574 dest.first, dest.second, norm);
2577 template <
class SrcTraverser,
class SrcAccessor,
2578 class DestTraverser,
class DestAccessor>
2580 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2581 DestTraverser dul, DestAccessor dest, fftw_real norm)
2583 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2584 norm, FFTW_RODFT00, FFTW_REDFT00);
2587 template <
class DestTraverser,
class DestAccessor>
2589 fourierTransformRealOE(
2593 DestTraverser dul, DestAccessor dest, fftw_real norm)
2595 int w = slr.x - sul.x;
2598 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2599 fourierTransformRealImpl(sul, slr, dul, dest,
2600 norm, FFTW_RODFT00, FFTW_REDFT00);
2602 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2603 norm, FFTW_RODFT00, FFTW_REDFT00);
2608 template <
class SrcTraverser,
class SrcAccessor,
2609 class DestTraverser,
class DestAccessor>
2611 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2612 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2614 fourierTransformRealEO(src.first, src.second, src.third,
2615 dest.first, dest.second, norm);
2618 template <
class SrcTraverser,
class SrcAccessor,
2619 class DestTraverser,
class DestAccessor>
2621 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2622 DestTraverser dul, DestAccessor dest, fftw_real norm)
2624 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2625 norm, FFTW_REDFT00, FFTW_RODFT00);
2628 template <
class DestTraverser,
class DestAccessor>
2630 fourierTransformRealEO(
2634 DestTraverser dul, DestAccessor dest, fftw_real norm)
2636 int w = slr.x - sul.x;
2639 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2640 fourierTransformRealImpl(sul, slr, dul, dest,
2641 norm, FFTW_REDFT00, FFTW_RODFT00);
2643 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2644 norm, FFTW_REDFT00, FFTW_RODFT00);
2649 template <
class SrcTraverser,
class SrcAccessor,
2650 class DestTraverser,
class DestAccessor>
2652 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2653 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2655 fourierTransformRealOO(src.first, src.second, src.third,
2656 dest.first, dest.second, norm);
2659 template <
class SrcTraverser,
class SrcAccessor,
2660 class DestTraverser,
class DestAccessor>
2662 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2663 DestTraverser dul, DestAccessor dest, fftw_real norm)
2665 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2666 norm, FFTW_RODFT00, FFTW_RODFT00);
2669 template <
class DestTraverser,
class DestAccessor>
2671 fourierTransformRealOO(
2675 DestTraverser dul, DestAccessor dest, fftw_real norm)
2677 int w = slr.x - sul.x;
2680 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2681 fourierTransformRealImpl(sul, slr, dul, dest,
2682 norm, FFTW_RODFT00, FFTW_RODFT00);
2684 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2685 norm, FFTW_RODFT00, FFTW_RODFT00);
2690 template <
class SrcTraverser,
class SrcAccessor,
2691 class DestTraverser,
class DestAccessor>
2693 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2694 DestTraverser dul, DestAccessor dest,
2695 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2698 copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2700 fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
2701 dul, dest,
norm, kindx, kindy);
2705 template <
class DestTraverser,
class DestAccessor>
2707 fourierTransformRealImpl(
2710 DestTraverser dul, DestAccessor dest,
2711 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2713 int w = slr.x - sul.x;
2714 int h = slr.y - sul.y;
2715 BasicImage<fftw_real> res(w, h);
2717 fftw_plan plan = fftw_plan_r2r_2d(h, w,
2718 (fftw_real *)&(*sul), (fftw_real *)res.begin(),
2719 kindy, kindx, FFTW_ESTIMATE);
2721 fftw_destroy_plan(plan);
2725 std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2727 copyImage(srcImageRange(res), destIter(dul, dest));
2735 #endif // VIGRA_FFTW3_HXX