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>
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>
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>
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>
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)
2088 template <
class SrcImageIterator,
class SrcAccessor,
2089 class FilterImageIterator,
class FilterAccessor,
2090 class DestImageIterator,
class DestAccessor>
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);
2101 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2102 destImage(workImage, FFTWWriteRealAccessor<>()));
2106 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2107 filterUpperLeft, fa,
2111 template <
class FilterImageIterator,
class FilterAccessor,
2112 class DestImageIterator,
class DestAccessor>
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,
2132 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2133 destImage(workImage));
2136 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2137 filterUpperLeft, fa,
2142 template <
class SrcImageIterator,
class SrcAccessor,
2143 class FilterImageIterator,
class FilterAccessor,
2144 class DestImageIterator,
class DestAccessor>
2147 pair<FilterImageIterator, FilterAccessor> filter,
2148 pair<DestImageIterator, DestAccessor> dest)
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);
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>
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>
2319 const ImageArray<FilterType> &filters,
2320 ImageArray<DestImage> &results)
2326 template <
class SrcImageIterator,
class SrcAccessor,
2327 class FilterType,
class DestImage>
2329 SrcImageIterator srcLowerRight, SrcAccessor sa,
2330 const ImageArray<FilterType> &filters,
2331 ImageArray<DestImage> &results)
2333 int w = int(srcLowerRight.x - srcUpperLeft.x);
2334 int h = int(srcLowerRight.y - srcUpperLeft.y);
2337 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2338 destImage(workImage, FFTWWriteRealAccessor<>()));
2341 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2345 template <
class FilterType,
class DestImage>
2351 const ImageArray<FilterType> &filters,
2352 ImageArray<DestImage> &results)
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;
2364 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2365 destImage(workImage));
2368 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2373 template <
class FilterType,
class DestImage>
2374 void applyFourierFilterFamilyImpl(
2378 const ImageArray<FilterType> &filters,
2379 ImageArray<DestImage> &results)
2385 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
2386 "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2389 results.resize(filters.size());
2390 results.resizeImages(filters.imageSize());
2393 int w = int(srcLowerRight.x - srcUpperLeft.x);
2394 int h = int(srcLowerRight.y - srcUpperLeft.y);
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++)
2418 destImage(result), std::multiplies<FFTWComplex<> >());
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)
2735 copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2737 fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
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;
2752 BasicImage<fftw_real> res(w, h);
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