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>& ) throw()
699 template<
class Other>
700 FFTWAllocator& operator=(
const FFTWAllocator<Other>& )
705 pointer allocate(size_type count,
void * = 0)
707 return (pointer)fftw_malloc(count *
sizeof(Ty));
710 void deallocate(pointer ptr, size_type )
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>& ) throw()
767 template<
class Other>
768 allocator& operator=(
const allocator<Other>& )
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 std::log((*i).magnitude() + 1);
1416 template <
class ITERATOR,
class DIFFERENCE>
1418 return std::log((i[d]).magnitude() + 1);
1428 template <
class Real =
double>
1436 template <
class ITERATOR>
1438 return (*i).phase();
1442 template <
class ITERATOR,
class DIFFERENCE>
1444 return (i[d]).phase();
1620 template <
class SrcImageIterator,
class SrcAccessor,
1621 class DestImageIterator,
class DestAccessor>
1623 SrcImageIterator src_lowerright, SrcAccessor sa,
1624 DestImageIterator dest_upperleft, DestAccessor da)
1626 int w = int(src_lowerright.x - src_upperleft.x);
1627 int h = int(src_lowerright.y - src_upperleft.y);
1635 src_upperleft + Diff2D(w2, h2), sa),
1636 destIter (dest_upperleft + Diff2D(w1, h1), da));
1639 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1640 src_lowerright, sa),
1641 destIter (dest_upperleft, da));
1644 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1645 src_upperleft + Diff2D(w, h2), sa),
1646 destIter (dest_upperleft + Diff2D(0, h1), da));
1649 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1650 src_upperleft + Diff2D(w2, h), sa),
1651 destIter (dest_upperleft + Diff2D(w1, 0), da));
1654 template <
class SrcImageIterator,
class SrcAccessor,
1655 class DestImageIterator,
class DestAccessor>
1657 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1658 pair<DestImageIterator, DestAccessor> dest)
1661 dest.first, dest.second);
1711 template <
class SrcImageIterator,
class SrcAccessor,
1712 class DestImageIterator,
class DestAccessor>
1714 SrcImageIterator src_lowerright, SrcAccessor sa,
1715 DestImageIterator dest_upperleft, DestAccessor da)
1717 int w = int(src_lowerright.x - src_upperleft.x);
1718 int h = int(src_lowerright.y - src_upperleft.y);
1726 src_upperleft + Diff2D(w2, h2), sa),
1727 destIter (dest_upperleft + Diff2D(w1, h1), da));
1730 copyImage(srcIterRange(src_upperleft + Diff2D(w2, h2),
1731 src_lowerright, sa),
1732 destIter (dest_upperleft, da));
1735 copyImage(srcIterRange(src_upperleft + Diff2D(w2, 0),
1736 src_upperleft + Diff2D(w, h2), sa),
1737 destIter (dest_upperleft + Diff2D(0, h1), da));
1740 copyImage(srcIterRange(src_upperleft + Diff2D(0, h2),
1741 src_upperleft + Diff2D(w2, h), sa),
1742 destIter (dest_upperleft + Diff2D(w1, 0), da));
1745 template <
class SrcImageIterator,
class SrcAccessor,
1746 class DestImageIterator,
class DestAccessor>
1748 triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1749 pair<DestImageIterator, DestAccessor> dest)
1752 dest.first, dest.second);
1763 int w = int(slr.x - sul.x);
1764 int h = int(slr.y - sul.y);
1768 fftw_complex * srcPtr = (fftw_complex *)(&*sul);
1769 fftw_complex * destPtr = (fftw_complex *)(&*dul);
1772 if (h > 1 && &(*(sul + Diff2D(w, 0))) != &(*(sul + Diff2D(0, 1))))
1774 sworkImage.resize(w, h);
1775 copyImage(srcIterRange(sul, slr, src), destImage(sworkImage));
1776 srcPtr = (fftw_complex *)(&(*sworkImage.upperLeft()));
1778 if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1780 dworkImage.resize(w, h);
1781 destPtr = (fftw_complex *)(&(*dworkImage.upperLeft()));
1784 fftw_plan plan = fftw_plan_dft_2d(h, w, srcPtr, destPtr, sign, FFTW_ESTIMATE );
1786 fftw_destroy_plan(plan);
1788 if (h > 1 && &(*(dul + Diff2D(w, 0))) != &(*(dul + Diff2D(0, 1))))
1790 copyImage(srcImageRange(dworkImage), destIter(dul, dest));
1949 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_FORWARD);
1952 template <
class SrcImageIterator,
class SrcAccessor>
1954 SrcImageIterator srcLowerRight, SrcAccessor sa,
1958 int w= srcLowerRight.x - srcUpperLeft.x;
1959 int h= srcLowerRight.y - srcUpperLeft.y;
1962 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
1963 destImage(workImage, FFTWWriteRealAccessor<>()));
1967 fourierTransform(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
1971 template <
class SrcImageIterator,
class SrcAccessor>
1973 void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
1974 pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest)
1976 fourierTransform(src.first, src.second, src.third, dest.first, dest.second);
1990 detail::fourierTransformImpl(sul, slr, src, dul, dest, FFTW_BACKWARD);
1993 template <
class DestImageIterator,
class DestAccessor>
1996 DestImageIterator dul, DestAccessor dest)
1998 int w = slr.x - sul.x;
1999 int h = slr.y - sul.y;
2003 copyImage(srcImageRange(workImage), destIter(dul, dest));
2007 template <
class DestImageIterator,
class DestAccessor>
2011 pair<DestImageIterator, DestAccessor> dest)
2114 template <
class SrcImageIterator,
class SrcAccessor,
2115 class FilterImageIterator,
class FilterAccessor,
2116 class DestImageIterator,
class DestAccessor>
2118 SrcImageIterator srcLowerRight, SrcAccessor sa,
2119 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2120 DestImageIterator destUpperLeft, DestAccessor da)
2123 int w = int(srcLowerRight.x - srcUpperLeft.x);
2124 int h = int(srcLowerRight.y - srcUpperLeft.y);
2127 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2128 destImage(workImage, FFTWWriteRealAccessor<>()));
2132 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2133 filterUpperLeft, fa,
2137 template <
class FilterImageIterator,
class FilterAccessor,
2138 class DestImageIterator,
class DestAccessor>
2144 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2145 DestImageIterator destUpperLeft, DestAccessor da)
2147 int w = srcLowerRight.x - srcUpperLeft.x;
2148 int h = srcLowerRight.y - srcUpperLeft.y;
2151 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2152 applyFourierFilterImpl(srcUpperLeft, srcLowerRight, sa,
2153 filterUpperLeft, fa,
2158 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2159 destImage(workImage));
2162 applyFourierFilterImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2163 filterUpperLeft, fa,
2168 template <
class SrcImageIterator,
class SrcAccessor,
2169 class FilterImageIterator,
class FilterAccessor,
2170 class DestImageIterator,
class DestAccessor>
2173 pair<FilterImageIterator, FilterAccessor> filter,
2174 pair<DestImageIterator, DestAccessor> dest)
2177 filter.first, filter.second,
2178 dest.first, dest.second);
2181 template <
class FilterImageIterator,
class FilterAccessor,
2182 class DestImageIterator,
class DestAccessor>
2183 void applyFourierFilterImpl(
2187 FilterImageIterator filterUpperLeft, FilterAccessor fa,
2188 DestImageIterator destUpperLeft, DestAccessor da)
2190 int w = int(srcLowerRight.x - srcUpperLeft.x);
2191 int h = int(srcLowerRight.y - srcUpperLeft.y);
2196 fftw_plan forwardPlan=
2197 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2198 (fftw_complex *)complexResultImg.begin(),
2199 FFTW_FORWARD, FFTW_ESTIMATE );
2200 fftw_execute(forwardPlan);
2201 fftw_destroy_plan(forwardPlan);
2204 combineTwoImages(srcImageRange(complexResultImg), srcIter(filterUpperLeft, fa),
2205 destImage(complexResultImg), std::multiplies<FFTWComplex<> >());
2208 fftw_plan backwardPlan=
2209 fftw_plan_dft_2d(h, w, (fftw_complex *)complexResultImg.begin(),
2210 (fftw_complex *)complexResultImg.begin(),
2211 FFTW_BACKWARD, FFTW_ESTIMATE);
2212 fftw_execute(backwardPlan);
2213 fftw_destroy_plan(backwardPlan);
2216 NumericTraits<typename DestAccessor::value_type>::isScalar
2220 applyFourierFilterImplNormalization(complexResultImg, destUpperLeft, da,
2224 template <
class DestImageIterator,
class DestAccessor>
2226 DestImageIterator destUpperLeft,
2230 double normFactor= 1.0/(srcImage.width() * srcImage.height());
2232 for(
int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2234 DestImageIterator dIt= destUpperLeft;
2235 for(
int x= 0; x< srcImage.width(); x++, dIt.x++)
2237 da.setComponent(srcImage(x, y).re()*normFactor, dIt, 0);
2238 da.setComponent(srcImage(x, y).im()*normFactor, dIt, 1);
2244 void applyFourierFilterImplNormalization(
FFTWComplexImage const & srcImage,
2249 transformImage(srcImageRange(srcImage), destIter(destUpperLeft, da),
2253 template <
class DestImageIterator,
class DestAccessor>
2254 void applyFourierFilterImplNormalization(
FFTWComplexImage const & srcImage,
2255 DestImageIterator destUpperLeft,
2259 double normFactor= 1.0/(srcImage.width() * srcImage.height());
2261 for(
int y=0; y<srcImage.height(); y++, destUpperLeft.y++)
2263 DestImageIterator dIt= destUpperLeft;
2264 for(
int x= 0; x< srcImage.width(); x++, dIt.x++)
2265 da.set(srcImage(x, y).re()*normFactor, dIt);
2341 template <
class SrcImageIterator,
class SrcAccessor,
2342 class FilterType,
class DestImage>
2345 const ImageArray<FilterType> &filters,
2346 ImageArray<DestImage> &results)
2352 template <
class SrcImageIterator,
class SrcAccessor,
2353 class FilterType,
class DestImage>
2355 SrcImageIterator srcLowerRight, SrcAccessor sa,
2356 const ImageArray<FilterType> &filters,
2357 ImageArray<DestImage> &results)
2359 int w = int(srcLowerRight.x - srcUpperLeft.x);
2360 int h = int(srcLowerRight.y - srcUpperLeft.y);
2363 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2364 destImage(workImage, FFTWWriteRealAccessor<>()));
2367 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2371 template <
class FilterType,
class DestImage>
2377 const ImageArray<FilterType> &filters,
2378 ImageArray<DestImage> &results)
2380 int w= srcLowerRight.x - srcUpperLeft.x;
2383 if (&(*(srcUpperLeft + Diff2D(w, 0))) == &(*(srcUpperLeft + Diff2D(0, 1))))
2384 applyFourierFilterFamilyImpl(srcUpperLeft, srcLowerRight, sa,
2388 int h = srcLowerRight.y - srcUpperLeft.y;
2390 copyImage(srcIterRange(srcUpperLeft, srcLowerRight, sa),
2391 destImage(workImage));
2394 applyFourierFilterFamilyImpl(cworkImage.upperLeft(), cworkImage.lowerRight(), cworkImage.accessor(),
2399 template <
class FilterType,
class DestImage>
2400 void applyFourierFilterFamilyImpl(
2404 const ImageArray<FilterType> &filters,
2405 ImageArray<DestImage> &results)
2411 vigra_precondition((srcLowerRight - srcUpperLeft) == filters.imageSize(),
2412 "applyFourierFilterFamily called with src image size != filters.imageSize()!");
2415 results.resize(filters.size());
2416 results.resizeImages(filters.imageSize());
2419 int w = int(srcLowerRight.x - srcUpperLeft.x);
2420 int h = int(srcLowerRight.y - srcUpperLeft.y);
2425 fftw_plan forwardPlan=
2426 fftw_plan_dft_2d(h, w, (fftw_complex *)&(*srcUpperLeft),
2427 (fftw_complex *)freqImage.begin(),
2428 FFTW_FORWARD, FFTW_ESTIMATE );
2429 fftw_execute(forwardPlan);
2430 fftw_destroy_plan(forwardPlan);
2432 fftw_plan backwardPlan=
2433 fftw_plan_dft_2d(h, w, (fftw_complex *)result.begin(),
2434 (fftw_complex *)result.begin(),
2435 FFTW_BACKWARD, FFTW_ESTIMATE );
2437 NumericTraits<typename DestImage::Accessor::value_type>::isScalar
2441 for (
unsigned int i= 0; i < filters.size(); i++)
2444 destImage(result), std::multiplies<FFTWComplex<> >());
2447 fftw_execute(backwardPlan);
2450 applyFourierFilterImplNormalization(result,
2451 results[i].upperLeft(), results[i].accessor(),
2454 fftw_destroy_plan(backwardPlan);
2589 template <
class SrcTraverser,
class SrcAccessor,
2590 class DestTraverser,
class DestAccessor>
2592 fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2593 pair<DestTraverser, DestAccessor> dest, fftw_real
norm)
2595 fourierTransformRealEE(src.first, src.second, src.third,
2596 dest.first, dest.second, norm);
2599 template <
class SrcTraverser,
class SrcAccessor,
2600 class DestTraverser,
class DestAccessor>
2602 fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2603 DestTraverser dul, DestAccessor dest, fftw_real
norm)
2605 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2606 norm, FFTW_REDFT00, FFTW_REDFT00);
2609 template <
class DestTraverser,
class DestAccessor>
2611 fourierTransformRealEE(
2615 DestTraverser dul, DestAccessor dest, fftw_real
norm)
2617 int w = slr.x - sul.x;
2620 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2621 fourierTransformRealImpl(sul, slr, dul, dest,
2622 norm, FFTW_REDFT00, FFTW_REDFT00);
2624 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2625 norm, FFTW_REDFT00, FFTW_REDFT00);
2630 template <
class SrcTraverser,
class SrcAccessor,
2631 class DestTraverser,
class DestAccessor>
2633 fourierTransformRealOE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2634 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2636 fourierTransformRealOE(src.first, src.second, src.third,
2637 dest.first, dest.second, norm);
2640 template <
class SrcTraverser,
class SrcAccessor,
2641 class DestTraverser,
class DestAccessor>
2643 fourierTransformRealOE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2644 DestTraverser dul, DestAccessor dest, fftw_real norm)
2646 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2647 norm, FFTW_RODFT00, FFTW_REDFT00);
2650 template <
class DestTraverser,
class DestAccessor>
2652 fourierTransformRealOE(
2656 DestTraverser dul, DestAccessor dest, fftw_real norm)
2658 int w = slr.x - sul.x;
2661 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2662 fourierTransformRealImpl(sul, slr, dul, dest,
2663 norm, FFTW_RODFT00, FFTW_REDFT00);
2665 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2666 norm, FFTW_RODFT00, FFTW_REDFT00);
2671 template <
class SrcTraverser,
class SrcAccessor,
2672 class DestTraverser,
class DestAccessor>
2674 fourierTransformRealEO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2675 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2677 fourierTransformRealEO(src.first, src.second, src.third,
2678 dest.first, dest.second, norm);
2681 template <
class SrcTraverser,
class SrcAccessor,
2682 class DestTraverser,
class DestAccessor>
2684 fourierTransformRealEO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2685 DestTraverser dul, DestAccessor dest, fftw_real norm)
2687 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2688 norm, FFTW_REDFT00, FFTW_RODFT00);
2691 template <
class DestTraverser,
class DestAccessor>
2693 fourierTransformRealEO(
2697 DestTraverser dul, DestAccessor dest, fftw_real norm)
2699 int w = slr.x - sul.x;
2702 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2703 fourierTransformRealImpl(sul, slr, dul, dest,
2704 norm, FFTW_REDFT00, FFTW_RODFT00);
2706 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2707 norm, FFTW_REDFT00, FFTW_RODFT00);
2712 template <
class SrcTraverser,
class SrcAccessor,
2713 class DestTraverser,
class DestAccessor>
2715 fourierTransformRealOO(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
2716 pair<DestTraverser, DestAccessor> dest, fftw_real norm)
2718 fourierTransformRealOO(src.first, src.second, src.third,
2719 dest.first, dest.second, norm);
2722 template <
class SrcTraverser,
class SrcAccessor,
2723 class DestTraverser,
class DestAccessor>
2725 fourierTransformRealOO(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2726 DestTraverser dul, DestAccessor dest, fftw_real norm)
2728 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2729 norm, FFTW_RODFT00, FFTW_RODFT00);
2732 template <
class DestTraverser,
class DestAccessor>
2734 fourierTransformRealOO(
2738 DestTraverser dul, DestAccessor dest, fftw_real norm)
2740 int w = slr.x - sul.x;
2743 if (&(*(sul + Diff2D(w, 0))) == &(*(sul + Diff2D(0, 1))))
2744 fourierTransformRealImpl(sul, slr, dul, dest,
2745 norm, FFTW_RODFT00, FFTW_RODFT00);
2747 fourierTransformRealWorkImageImpl(sul, slr, src, dul, dest,
2748 norm, FFTW_RODFT00, FFTW_RODFT00);
2753 template <
class SrcTraverser,
class SrcAccessor,
2754 class DestTraverser,
class DestAccessor>
2756 fourierTransformRealWorkImageImpl(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
2757 DestTraverser dul, DestAccessor dest,
2758 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2761 copyImage(srcIterRange(sul, slr, src), destImage(workImage));
2763 fourierTransformRealImpl(cworkImage.upperLeft(), cworkImage.lowerRight(),
2764 dul, dest,
norm, kindx, kindy);
2768 template <
class DestTraverser,
class DestAccessor>
2770 fourierTransformRealImpl(
2773 DestTraverser dul, DestAccessor dest,
2774 fftw_real norm, fftw_r2r_kind kindx, fftw_r2r_kind kindy)
2776 int w = slr.x - sul.x;
2777 int h = slr.y - sul.y;
2778 BasicImage<fftw_real> res(w, h);
2780 fftw_plan plan = fftw_plan_r2r_2d(h, w,
2781 (fftw_real *)&(*sul), (fftw_real *)res.begin(),
2782 kindy, kindx, FFTW_ESTIMATE);
2784 fftw_destroy_plan(plan);
2788 std::bind1st(std::multiplies<fftw_real>(), 1.0 / norm));
2790 copyImage(srcImageRange(res), destIter(dul, dest));
2798 #endif // VIGRA_FFTW3_HXX
FFTWComplex(fftwf_complex const &o)
Definition: fftw3.hxx:202
Definition: fftw3.hxx:1217
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Definition: fftw3.hxx:1339
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
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:1443
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
void set(FFTWComplex< R > const &v, ITERATOR const &i) const
Write imaginary part at iterator position into a scalar.
Definition: fftw3.hxx:1300
void set(value_type const &v, ITERATOR const &i) const
Definition: fftw3.hxx:1330
value_type SquaredNormType
Definition: fftw3.hxx:164
R arg(const FFTWComplex< R > &a)
phase
Definition: fftw3.hxx:1009
FFTWReal2Complex< Real >::type complex_type
Definition: fftw3.hxx:136
Definition: fftw3.hxx:1351
SquaredNormType squaredMagnitude() const
Definition: fftw3.hxx:358
Diff2D operator-(Diff2D const &a, Diff2D const &b)
Definition: diff2d.hxx:711
void set(FFTWComplex< R > const &v, ITERATOR const &i) const
Write real part at iterator position into a scalar.
Definition: fftw3.hxx:1250
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
IteratorTraits< traverser >::DefaultAccessor Accessor
Definition: basicimage.hxx:573
value_type phase() const
Definition: fftw3.hxx:368
void set(FFTWComplex< R > const &v, ITERATOR const &i, DIFFERENCE d) const
Write real part at offset from iterator position into a scalar.
Definition: fftw3.hxx:1256
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
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1355
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read natural log of magnitude at offset from iterator position.
Definition: fftw3.hxx:1417
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
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Write imaginary part at offset from iterator position from a scalar.
Definition: fftw3.hxx:1294
Diff2D operator+(Diff2D const &a, Diff2D const &b)
Definition: diff2d.hxx:739
void set(FFTWComplex< R > const &v, ITERATOR const &i, DIFFERENCE d) const
Write imaginary part at offset from iterator position into a scalar.
Definition: fftw3.hxx:1306
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
value_type operator()(ITERATOR const &i) const
Read squared magnitude at iterator position.
Definition: fftw3.hxx:1359
FFTWComplex< R > & operator-=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
subtract-assignment
Definition: fftw3.hxx:867
BasicImageIterator< PIXELTYPE, PIXELTYPE ** > traverser
Definition: basicimage.hxx:528
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read squared magnitude at offset from iterator position.
Definition: fftw3.hxx:1365
Definition: fftw3.hxx:1429
value_type operator()(ITERATOR const &i) const
Read natural log of magnitude at iterator position.
Definition: fftw3.hxx:1411
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 & operator=(float o)
Definition: fftw3.hxx:291
Definition: basicimage.hxx:262
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read real part at offset from iterator position.
Definition: fftw3.hxx:1232
FFTWComplex< R > & operator+=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
add-assignment
Definition: fftw3.hxx:859
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
bool operator!=(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
not equal
Definition: fftw3.hxx:841
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:1437
FFTWComplex & operator=(long double o)
Definition: fftw3.hxx:300
bool operator==(FFTWComplex< R > const &a, const FFTWComplex< R > &b)
equal
Definition: fftw3.hxx:825
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
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void copyImage(...)
Copy source image into destination image.
FFTWComplex< R > & operator*=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
multiply-assignment
Definition: fftw3.hxx:875
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
Class for fixed size vectors.This class contains an array of size SIZE of the specified VALUETYPE...
Definition: accessor.hxx:940
IteratorTraits< const_traverser >::DefaultAccessor ConstAccessor
Definition: basicimage.hxx:578
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
ConstBasicImageIterator< PIXELTYPE, PIXELTYPE ** > const_traverser
Definition: basicimage.hxx:538
Definition: fftw3.hxx:1377
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
Fundamental class template for images.
Definition: basicimage.hxx:475
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1407
void set(value_type const &v, ITERATOR const &i) const
Write imaginary part at iterator position from a scalar.
Definition: fftw3.hxx:1288
FFTWComplex operator-() const
Definition: fftw3.hxx:353
Definition: basicimage.hxx:294
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
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
FFTWComplex< R > & operator/=(FFTWComplex< R > &a, const FFTWComplex< R > &b)
divide-assignment
Definition: fftw3.hxx:884
NormType magnitude() const
Definition: fftw3.hxx:363
void set(value_type const &v, ITERATOR const &i) const
Write real part at iterator position from a scalar.
Definition: fftw3.hxx:1238
Definition: fftw3.hxx:1403
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
T sign(T t)
The sign function.
Definition: mathutil.hxx:591
Real value_type
The accessor's value type.
Definition: fftw3.hxx:1433
value_type operator()(ITERATOR const &i, DIFFERENCE d) const
Read magnitude at offset from iterator position.
Definition: fftw3.hxx:1391
void set(value_type const &v, ITERATOR const &i, DIFFERENCE d) const
Write real part at offset from iterator position from a scalar.
Definition: fftw3.hxx:1244
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(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
FFTWComplex & operator=(fftwf_complex const &o)
Definition: fftw3.hxx:264