37 #ifndef VIGRA_SEPARABLECONVOLUTION_HXX
38 #define VIGRA_SEPARABLECONVOLUTION_HXX
41 #include "utilities.hxx"
42 #include "numerictraits.hxx"
43 #include "imageiteratoradapter.hxx"
44 #include "bordertreatment.hxx"
45 #include "gaussians.hxx"
46 #include "array_vector.hxx"
56 template <
class SrcIterator,
class SrcAccessor,
57 class DestIterator,
class DestAccessor,
58 class KernelIterator,
class KernelAccessor>
59 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
60 DestIterator
id, DestAccessor da,
61 KernelIterator kernel, KernelAccessor ka,
62 int kleft,
int kright)
65 int w = std::distance( is, iend );
67 typedef typename PromoteTraits<
68 typename SrcAccessor::value_type,
69 typename KernelAccessor::value_type>::Promote SumType;
71 SrcIterator ibegin = is;
73 for(
int x=0; x<w; ++x, ++is, ++id)
75 KernelIterator ik = kernel + kright;
76 SumType sum = NumericTraits<SumType>::zero();
81 SrcIterator iss = iend + x0;
83 for(; x0; ++x0, --ik, ++iss)
85 sum += ka(ik) * sa(iss);
89 SrcIterator isend = is + (1 - kleft);
90 for(; iss != isend ; --ik, ++iss)
92 sum += ka(ik) * sa(iss);
95 else if(w-x <= -kleft)
97 SrcIterator iss = is + (-kright);
98 SrcIterator isend = iend;
99 for(; iss != isend ; --ik, ++iss)
101 sum += ka(ik) * sa(iss);
104 int x0 = -kleft - w + x + 1;
107 for(; x0; --x0, --ik, ++iss)
109 sum += ka(ik) * sa(iss);
114 SrcIterator iss = is - kright;
115 SrcIterator isend = is + (1 - kleft);
116 for(; iss != isend ; --ik, ++iss)
118 sum += ka(ik) * sa(iss);
122 da.set(detail::RequiresExplicitCast<
typename
123 DestAccessor::value_type>::cast(sum),
id);
133 template <
class SrcIterator,
class SrcAccessor,
134 class DestIterator,
class DestAccessor,
135 class KernelIterator,
class KernelAccessor,
137 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
138 DestIterator
id, DestAccessor da,
139 KernelIterator kernel, KernelAccessor ka,
140 int kleft,
int kright, Norm
norm)
143 int w = std::distance( is, iend );
145 typedef typename PromoteTraits<
146 typename SrcAccessor::value_type,
147 typename KernelAccessor::value_type>::Promote SumType;
149 SrcIterator ibegin = is;
151 for(
int x=0; x<w; ++x, ++is, ++id)
153 KernelIterator ik = kernel + kright;
154 SumType sum = NumericTraits<SumType>::zero();
159 Norm clipped = NumericTraits<Norm>::zero();
161 for(; x0; ++x0, --ik)
166 SrcIterator iss = ibegin;
167 SrcIterator isend = is + (1 - kleft);
168 for(; iss != isend ; --ik, ++iss)
170 sum += ka(ik) * sa(iss);
173 sum = norm / (norm - clipped) * sum;
175 else if(w-x <= -kleft)
177 SrcIterator iss = is + (-kright);
178 SrcIterator isend = iend;
179 for(; iss != isend ; --ik, ++iss)
181 sum += ka(ik) * sa(iss);
184 Norm clipped = NumericTraits<Norm>::zero();
186 int x0 = -kleft - w + x + 1;
188 for(; x0; --x0, --ik)
193 sum = norm / (norm - clipped) * sum;
197 SrcIterator iss = is + (-kright);
198 SrcIterator isend = is + (1 - kleft);
199 for(; iss != isend ; --ik, ++iss)
201 sum += ka(ik) * sa(iss);
205 da.set(detail::RequiresExplicitCast<
typename
206 DestAccessor::value_type>::cast(sum),
id);
216 template <
class SrcIterator,
class SrcAccessor,
217 class DestIterator,
class DestAccessor,
218 class KernelIterator,
class KernelAccessor>
219 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
220 DestIterator
id, DestAccessor da,
221 KernelIterator kernel, KernelAccessor ka,
222 int kleft,
int kright)
225 int w = std::distance( is, iend );
227 typedef typename PromoteTraits<
228 typename SrcAccessor::value_type,
229 typename KernelAccessor::value_type>::Promote SumType;
231 SrcIterator ibegin = is;
233 for(
int x=0; x<w; ++x, ++is, ++id)
235 KernelIterator ik = kernel + kright;
236 SumType sum = NumericTraits<SumType>::zero();
241 SrcIterator iss = ibegin - x0;
243 for(; x0; ++x0, --ik, --iss)
245 sum += ka(ik) * sa(iss);
248 SrcIterator isend = is + (1 - kleft);
249 for(; iss != isend ; --ik, ++iss)
251 sum += ka(ik) * sa(iss);
254 else if(w-x <= -kleft)
256 SrcIterator iss = is + (-kright);
257 SrcIterator isend = iend;
258 for(; iss != isend ; --ik, ++iss)
260 sum += ka(ik) * sa(iss);
263 int x0 = -kleft - w + x + 1;
266 for(; x0; --x0, --ik, --iss)
268 sum += ka(ik) * sa(iss);
273 SrcIterator iss = is + (-kright);
274 SrcIterator isend = is + (1 - kleft);
275 for(; iss != isend ; --ik, ++iss)
277 sum += ka(ik) * sa(iss);
281 da.set(detail::RequiresExplicitCast<
typename
282 DestAccessor::value_type>::cast(sum),
id);
292 template <
class SrcIterator,
class SrcAccessor,
293 class DestIterator,
class DestAccessor,
294 class KernelIterator,
class KernelAccessor>
295 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
296 DestIterator
id, DestAccessor da,
297 KernelIterator kernel, KernelAccessor ka,
298 int kleft,
int kright)
301 int w = std::distance( is, iend );
303 typedef typename PromoteTraits<
304 typename SrcAccessor::value_type,
305 typename KernelAccessor::value_type>::Promote SumType;
307 SrcIterator ibegin = is;
309 for(
int x=0; x<w; ++x, ++is, ++id)
311 KernelIterator ik = kernel + kright;
312 SumType sum = NumericTraits<SumType>::zero();
317 SrcIterator iss = ibegin;
319 for(; x0; ++x0, --ik)
321 sum += ka(ik) * sa(iss);
324 SrcIterator isend = is + (1 - kleft);
325 for(; iss != isend ; --ik, ++iss)
327 sum += ka(ik) * sa(iss);
330 else if(w-x <= -kleft)
332 SrcIterator iss = is + (-kright);
333 SrcIterator isend = iend;
334 for(; iss != isend ; --ik, ++iss)
336 sum += ka(ik) * sa(iss);
339 int x0 = -kleft - w + x + 1;
342 for(; x0; --x0, --ik)
344 sum += ka(ik) * sa(iss);
349 SrcIterator iss = is + (-kright);
350 SrcIterator isend = is + (1 - kleft);
351 for(; iss != isend ; --ik, ++iss)
353 sum += ka(ik) * sa(iss);
357 da.set(detail::RequiresExplicitCast<
typename
358 DestAccessor::value_type>::cast(sum),
id);
368 template <
class SrcIterator,
class SrcAccessor,
369 class DestIterator,
class DestAccessor,
370 class KernelIterator,
class KernelAccessor>
371 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
372 DestIterator
id, DestAccessor da,
373 KernelIterator kernel, KernelAccessor ka,
374 int kleft,
int kright)
377 int w = std::distance( is, iend );
379 typedef typename PromoteTraits<
380 typename SrcAccessor::value_type,
381 typename KernelAccessor::value_type>::Promote SumType;
386 for(
int x=kright; x<w+kleft; ++x, ++is, ++id)
388 KernelIterator ik = kernel + kright;
389 SumType sum = NumericTraits<SumType>::zero();
391 SrcIterator iss = is + (-kright);
392 SrcIterator isend = is + (1 - kleft);
393 for(; iss != isend ; --ik, ++iss)
395 sum += ka(ik) * sa(iss);
398 da.set(detail::RequiresExplicitCast<
typename
399 DestAccessor::value_type>::cast(sum),
id);
537 template <
class SrcIterator,
class SrcAccessor,
538 class DestIterator,
class DestAccessor,
539 class KernelIterator,
class KernelAccessor>
540 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
541 DestIterator
id, DestAccessor da,
542 KernelIterator ik, KernelAccessor ka,
543 int kleft,
int kright, BorderTreatmentMode border)
545 typedef typename KernelAccessor::value_type KernelValue;
547 vigra_precondition(kleft <= 0,
548 "convolveLine(): kleft must be <= 0.\n");
549 vigra_precondition(kright >= 0,
550 "convolveLine(): kright must be >= 0.\n");
553 int w = std::distance( is, iend );
555 vigra_precondition(w >= kright - kleft + 1,
556 "convolveLine(): kernel longer than line\n");
560 case BORDER_TREATMENT_WRAP:
562 internalConvolveLineWrap(is, iend, sa,
id, da, ik, ka, kleft, kright);
565 case BORDER_TREATMENT_AVOID:
567 internalConvolveLineAvoid(is, iend, sa,
id, da, ik, ka, kleft, kright);
570 case BORDER_TREATMENT_REFLECT:
572 internalConvolveLineReflect(is, iend, sa,
id, da, ik, ka, kleft, kright);
575 case BORDER_TREATMENT_REPEAT:
577 internalConvolveLineRepeat(is, iend, sa,
id, da, ik, ka, kleft, kright);
580 case BORDER_TREATMENT_CLIP:
583 typedef typename KernelAccessor::value_type KT;
584 KT norm = NumericTraits<KT>::zero();
585 KernelIterator iik = ik + kleft;
586 for(
int i=kleft; i<=kright; ++i, ++iik) norm += ka(iik);
588 vigra_precondition(norm != NumericTraits<KT>::zero(),
589 "convolveLine(): Norm of kernel must be != 0"
590 " in mode BORDER_TREATMENT_CLIP.\n");
592 internalConvolveLineClip(is, iend, sa,
id, da, ik, ka, kleft, kright, norm);
597 vigra_precondition(0,
598 "convolveLine(): Unknown border treatment mode.\n");
603 template <
class SrcIterator,
class SrcAccessor,
604 class DestIterator,
class DestAccessor,
605 class KernelIterator,
class KernelAccessor>
607 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
608 pair<DestIterator, DestAccessor> dest,
609 tuple5<KernelIterator, KernelAccessor,
610 int,
int, BorderTreatmentMode> kernel)
613 dest.first, dest.second,
614 kernel.first, kernel.second,
615 kernel.third, kernel.fourth, kernel.fifth);
679 template <
class SrcIterator,
class SrcAccessor,
680 class DestIterator,
class DestAccessor,
681 class KernelIterator,
class KernelAccessor>
683 SrcIterator slowerright, SrcAccessor sa,
684 DestIterator dupperleft, DestAccessor da,
685 KernelIterator ik, KernelAccessor ka,
686 int kleft,
int kright, BorderTreatmentMode border)
688 typedef typename KernelAccessor::value_type KernelValue;
690 vigra_precondition(kleft <= 0,
691 "separableConvolveX(): kleft must be <= 0.\n");
692 vigra_precondition(kright >= 0,
693 "separableConvolveX(): kright must be >= 0.\n");
695 int w = slowerright.x - supperleft.x;
696 int h = slowerright.y - supperleft.y;
698 vigra_precondition(w >= kright - kleft + 1,
699 "separableConvolveX(): kernel longer than line\n");
703 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
705 typename SrcIterator::row_iterator rs = supperleft.rowIterator();
706 typename DestIterator::row_iterator rd = dupperleft.rowIterator();
709 ik, ka, kleft, kright, border);
713 template <
class SrcIterator,
class SrcAccessor,
714 class DestIterator,
class DestAccessor,
715 class KernelIterator,
class KernelAccessor>
718 pair<DestIterator, DestAccessor> dest,
719 tuple5<KernelIterator, KernelAccessor,
720 int,
int, BorderTreatmentMode> kernel)
723 dest.first, dest.second,
724 kernel.first, kernel.second,
725 kernel.third, kernel.fourth, kernel.fifth);
791 template <
class SrcIterator,
class SrcAccessor,
792 class DestIterator,
class DestAccessor,
793 class KernelIterator,
class KernelAccessor>
795 SrcIterator slowerright, SrcAccessor sa,
796 DestIterator dupperleft, DestAccessor da,
797 KernelIterator ik, KernelAccessor ka,
798 int kleft,
int kright, BorderTreatmentMode border)
800 typedef typename KernelAccessor::value_type KernelValue;
802 vigra_precondition(kleft <= 0,
803 "separableConvolveY(): kleft must be <= 0.\n");
804 vigra_precondition(kright >= 0,
805 "separableConvolveY(): kright must be >= 0.\n");
807 int w = slowerright.x - supperleft.x;
808 int h = slowerright.y - supperleft.y;
810 vigra_precondition(h >= kright - kleft + 1,
811 "separableConvolveY(): kernel longer than line\n");
815 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
817 typename SrcIterator::column_iterator cs = supperleft.columnIterator();
818 typename DestIterator::column_iterator cd = dupperleft.columnIterator();
821 ik, ka, kleft, kright, border);
825 template <
class SrcIterator,
class SrcAccessor,
826 class DestIterator,
class DestAccessor,
827 class KernelIterator,
class KernelAccessor>
830 pair<DestIterator, DestAccessor> dest,
831 tuple5<KernelIterator, KernelAccessor,
832 int,
int, BorderTreatmentMode> kernel)
835 dest.first, dest.second,
836 kernel.first, kernel.second,
837 kernel.third, kernel.fourth, kernel.fifth);
893 template <
class ARITHTYPE>
913 typedef typename InternalVector::iterator
Iterator;
917 typedef typename InternalVector::iterator
iterator;
934 : iter_(i), base_(i),
935 count_(count), sum_(count),
941 vigra_precondition(count_ == 1 || count_ == sum_,
942 "Kernel1D::initExplicitly(): "
943 "Wrong number of init values.");
969 static value_type one() {
return NumericTraits<value_type>::one(); }
979 border_treatment_(BORDER_TREATMENT_REFLECT),
982 kernel_.push_back(norm_);
988 : kernel_(k.kernel_),
991 border_treatment_(k.border_treatment_),
1014 border_treatment_ = k.border_treatment_;
1016 kernel_ = k.kernel_;
1039 int size = right_ - left_ + 1;
1040 for(
unsigned int i=0; i<kernel_.
size(); ++i) kernel_[i] = v;
1041 norm_ = (double)size*v;
1043 return InitProxy(kernel_.
begin(),
size, norm_);
1253 this->
initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1281 this->
initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1309 this->
initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1344 vigra_precondition(a >= 0.0 && a <= 0.125,
1345 "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
1535 this->
initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
1575 vigra_precondition(left <= 0,
1576 "Kernel1D::initExplicitly(): left border must be <= 0.");
1577 vigra_precondition(right >= 0,
1578 "Kernel1D::initExplicitly(): right border must be >= 0.");
1583 kernel_.resize(right - left + 1);
1616 return kernel_[location -
left()];
1621 return kernel_[location -
left()];
1634 int size()
const {
return right_ - left_ + 1; }
1639 {
return border_treatment_; }
1644 { border_treatment_ = new_mode; }
1673 InternalVector kernel_;
1675 BorderTreatmentMode border_treatment_;
1679 template <
class ARITHTYPE>
1681 unsigned int derivativeOrder,
1684 typedef typename NumericTraits<value_type>::RealPromote TmpType;
1688 TmpType sum = NumericTraits<TmpType>::zero();
1690 if(derivativeOrder == 0)
1692 for(; k < kernel_.end(); ++k)
1699 unsigned int faculty = 1;
1700 for(
unsigned int i = 2; i <= derivativeOrder; ++i)
1702 for(
double x = left() + offset; k < kernel_.end(); ++x, ++k)
1704 sum = TmpType(sum + *k *
VIGRA_CSTD::pow(-x,
int(derivativeOrder)) / faculty);
1708 vigra_precondition(sum != NumericTraits<value_type>::zero(),
1709 "Kernel1D<ARITHTYPE>::normalize(): "
1710 "Cannot normalize a kernel with sum = 0");
1713 k = kernel_.begin();
1714 for(; k != kernel_.end(); ++k)
1724 template <
class ARITHTYPE>
1728 vigra_precondition(std_dev >= 0.0,
1729 "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
1736 int radius = (int)(3.0 * std_dev + 0.5);
1741 kernel_.erase(kernel_.begin(), kernel_.end());
1742 kernel_.reserve(radius*2+1);
1744 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
1746 kernel_.push_back(gauss(x));
1753 kernel_.erase(kernel_.begin(), kernel_.end());
1754 kernel_.push_back(1.0);
1765 border_treatment_ = BORDER_TREATMENT_REFLECT;
1770 template <
class ARITHTYPE>
1774 vigra_precondition(std_dev >= 0.0,
1775 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
1780 int radius = (int)(3.0*std_dev + 0.5);
1784 double f = 2.0 / std_dev / std_dev;
1787 int maxIndex = (int)(2.0 * (radius + 5.0 *
VIGRA_CSTD::sqrt((
double)radius)) + 0.5);
1789 warray[maxIndex] = 0.0;
1790 warray[maxIndex-1] = 1.0;
1792 for(
int i = maxIndex-2; i >= radius; --i)
1794 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
1795 if(warray[i] > 1.0e40)
1797 warray[i+1] /= warray[i];
1805 warray[radius+1] = er * warray[radius+1] / warray[radius];
1806 warray[radius] = er;
1808 for(
int i = radius-1; i >= 0; --i)
1810 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
1814 double scale = norm / (2*er - warray[0]);
1816 initExplicitly(-radius, radius);
1819 for(
int i=0; i<=radius; ++i)
1821 c[i] = c[-i] = warray[i] * scale;
1826 kernel_.erase(kernel_.begin(), kernel_.end());
1827 kernel_.push_back(norm);
1835 border_treatment_ = BORDER_TREATMENT_REFLECT;
1840 template <
class ARITHTYPE>
1846 vigra_precondition(order >= 0,
1847 "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
1851 initGaussian(std_dev, norm);
1855 vigra_precondition(std_dev > 0.0,
1856 "Kernel1D::initGaussianDerivative(): "
1857 "Standard deviation must be > 0.");
1862 int radius = (int)(3.0 * std_dev + 0.5 * order + 0.5);
1868 kernel_.reserve(radius*2+1);
1873 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
1875 kernel_.push_back(gauss(x));
1876 dc += kernel_[kernel_.size()-1];
1878 dc = ARITHTYPE(dc / (2.0*radius + 1.0));
1884 for(
unsigned int i=0; i < kernel_.size(); ++i)
1894 normalize(norm, order);
1900 border_treatment_ = BORDER_TREATMENT_REFLECT;
1905 template <
class ARITHTYPE>
1910 vigra_precondition(radius > 0,
1911 "Kernel1D::initBinomial(): Radius must be > 0.");
1915 typename InternalVector::iterator x = kernel_.begin() + radius;
1919 for(
int j=radius-1; j>=-radius; --j)
1921 x[j] = 0.5 * x[j+1];
1922 for(
int i=j+1; i<radius; ++i)
1924 x[i] = 0.5 * (x[i] + x[i+1]);
1934 border_treatment_ = BORDER_TREATMENT_REFLECT;
1939 template <
class ARITHTYPE>
1943 vigra_precondition(radius > 0,
1944 "Kernel1D::initAveraging(): Radius must be > 0.");
1947 double scale = 1.0 / (radius * 2 + 1);
1950 kernel_.erase(kernel_.begin(), kernel_.end());
1951 kernel_.reserve(radius*2+1);
1953 for(
int i=0; i<=radius*2+1; ++i)
1955 kernel_.push_back(scale * norm);
1963 border_treatment_ = BORDER_TREATMENT_CLIP;
1968 template <
class ARITHTYPE>
1972 kernel_.erase(kernel_.begin(), kernel_.end());
1975 kernel_.push_back(ARITHTYPE(0.5 * norm));
1976 kernel_.push_back(ARITHTYPE(0.0 * norm));
1977 kernel_.push_back(ARITHTYPE(-0.5 * norm));
1985 border_treatment_ = BORDER_TREATMENT_REFLECT;
1996 template <
class KernelIterator,
class KernelAccessor>
1998 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
1999 kernel1d(KernelIterator ik, KernelAccessor ka,
2000 int kleft,
int kright, BorderTreatmentMode border)
2003 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2004 ik, ka, kleft, kright, border);
2009 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2010 int, int, BorderTreatmentMode>
2011 kernel1d(Kernel1D<T>
const & k)
2015 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2016 int, int, BorderTreatmentMode>(
2019 k.left(), k.right(),
2020 k.borderTreatment());
2025 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2026 int, int, BorderTreatmentMode>
2027 kernel1d(Kernel1D<T>
const & k, BorderTreatmentMode border)
2031 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2032 int, int, BorderTreatmentMode>(
2035 k.left(), k.right(),
2042 #endif // VIGRA_SEPARABLECONVOLUTION_HXX