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,
63 int start = 0,
int stop = 0)
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;
77 for(
int x=start; x<stop; ++x, ++is, ++id)
79 KernelIterator ik = kernel + kright;
80 SumType
sum = NumericTraits<SumType>::zero();
85 SrcIterator iss = iend + x0;
87 for(; x0; ++x0, --ik, ++iss)
89 sum += ka(ik) * sa(iss);
93 SrcIterator isend = is + (1 - kleft);
94 for(; iss != isend ; --ik, ++iss)
96 sum += ka(ik) * sa(iss);
99 else if(w-x <= -kleft)
101 SrcIterator iss = is + (-kright);
102 SrcIterator isend = iend;
103 for(; iss != isend ; --ik, ++iss)
105 sum += ka(ik) * sa(iss);
108 int x0 = -kleft - w + x + 1;
111 for(; x0; --x0, --ik, ++iss)
113 sum += ka(ik) * sa(iss);
118 SrcIterator iss = is - kright;
119 SrcIterator isend = is + (1 - kleft);
120 for(; iss != isend ; --ik, ++iss)
122 sum += ka(ik) * sa(iss);
126 da.set(detail::RequiresExplicitCast<
typename
127 DestAccessor::value_type>::cast(sum),
id);
137 template <
class SrcIterator,
class SrcAccessor,
138 class DestIterator,
class DestAccessor,
139 class KernelIterator,
class KernelAccessor,
141 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
142 DestIterator
id, DestAccessor da,
143 KernelIterator kernel, KernelAccessor ka,
144 int kleft,
int kright, Norm
norm,
145 int start = 0,
int stop = 0)
147 int w = std::distance( is, iend );
149 typedef typename PromoteTraits<
150 typename SrcAccessor::value_type,
151 typename KernelAccessor::value_type>::Promote SumType;
153 SrcIterator ibegin = is;
159 for(
int x=start; x<stop; ++x, ++is, ++id)
161 KernelIterator ik = kernel + kright;
162 SumType sum = NumericTraits<SumType>::zero();
167 Norm clipped = NumericTraits<Norm>::zero();
169 for(; x0; ++x0, --ik)
174 SrcIterator iss = ibegin;
175 SrcIterator isend = is + (1 - kleft);
176 for(; iss != isend ; --ik, ++iss)
178 sum += ka(ik) * sa(iss);
181 sum = norm / (norm - clipped) * sum;
183 else if(w-x <= -kleft)
185 SrcIterator iss = is + (-kright);
186 SrcIterator isend = iend;
187 for(; iss != isend ; --ik, ++iss)
189 sum += ka(ik) * sa(iss);
192 Norm clipped = NumericTraits<Norm>::zero();
194 int x0 = -kleft - w + x + 1;
196 for(; x0; --x0, --ik)
201 sum = norm / (norm - clipped) * sum;
205 SrcIterator iss = is + (-kright);
206 SrcIterator isend = is + (1 - kleft);
207 for(; iss != isend ; --ik, ++iss)
209 sum += ka(ik) * sa(iss);
213 da.set(detail::RequiresExplicitCast<
typename
214 DestAccessor::value_type>::cast(sum),
id);
224 template <
class SrcIterator,
class SrcAccessor,
225 class DestIterator,
class DestAccessor,
226 class KernelIterator,
class KernelAccessor>
227 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
228 DestIterator
id, DestAccessor da,
229 KernelIterator kernel, KernelAccessor ka,
230 int kleft,
int kright,
231 int start = 0,
int stop = 0)
233 int w = std::distance( is, iend );
235 typedef typename PromoteTraits<
236 typename SrcAccessor::value_type,
237 typename KernelAccessor::value_type>::Promote SumType;
239 SrcIterator ibegin = is;
245 for(
int x=start; x<stop; ++x, ++is, ++id)
247 KernelIterator ik = kernel + kright;
248 SumType sum = NumericTraits<SumType>::zero();
253 SrcIterator iss = ibegin - x0;
255 for(; x0; ++x0, --ik, --iss)
257 sum += ka(ik) * sa(iss);
260 SrcIterator isend = is + (1 - kleft);
261 for(; iss != isend ; --ik, ++iss)
263 sum += ka(ik) * sa(iss);
266 else if(w-x <= -kleft)
268 SrcIterator iss = is + (-kright);
269 SrcIterator isend = iend;
270 for(; iss != isend ; --ik, ++iss)
272 sum += ka(ik) * sa(iss);
275 int x0 = -kleft - w + x + 1;
278 for(; x0; --x0, --ik, --iss)
280 sum += ka(ik) * sa(iss);
285 SrcIterator iss = is + (-kright);
286 SrcIterator isend = is + (1 - kleft);
287 for(; iss != isend ; --ik, ++iss)
289 sum += ka(ik) * sa(iss);
293 da.set(detail::RequiresExplicitCast<
typename
294 DestAccessor::value_type>::cast(sum),
id);
304 template <
class SrcIterator,
class SrcAccessor,
305 class DestIterator,
class DestAccessor,
306 class KernelIterator,
class KernelAccessor>
307 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
308 DestIterator
id, DestAccessor da,
309 KernelIterator kernel, KernelAccessor ka,
310 int kleft,
int kright,
311 int start = 0,
int stop = 0)
313 int w = std::distance( is, iend );
315 typedef typename PromoteTraits<
316 typename SrcAccessor::value_type,
317 typename KernelAccessor::value_type>::Promote SumType;
319 SrcIterator ibegin = is;
325 for(
int x=start; x<stop; ++x, ++is, ++id)
327 KernelIterator ik = kernel + kright;
328 SumType sum = NumericTraits<SumType>::zero();
333 SrcIterator iss = ibegin;
335 for(; x0; ++x0, --ik)
337 sum += ka(ik) * sa(iss);
340 SrcIterator isend = is + (1 - kleft);
341 for(; iss != isend ; --ik, ++iss)
343 sum += ka(ik) * sa(iss);
346 else if(w-x <= -kleft)
348 SrcIterator iss = is + (-kright);
349 SrcIterator isend = iend;
350 for(; iss != isend ; --ik, ++iss)
352 sum += ka(ik) * sa(iss);
355 int x0 = -kleft - w + x + 1;
358 for(; x0; --x0, --ik)
360 sum += ka(ik) * sa(iss);
365 SrcIterator iss = is + (-kright);
366 SrcIterator isend = is + (1 - kleft);
367 for(; iss != isend ; --ik, ++iss)
369 sum += ka(ik) * sa(iss);
373 da.set(detail::RequiresExplicitCast<
typename
374 DestAccessor::value_type>::cast(sum),
id);
384 template <
class SrcIterator,
class SrcAccessor,
385 class DestIterator,
class DestAccessor,
386 class KernelIterator,
class KernelAccessor>
387 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
388 DestIterator
id, DestAccessor da,
389 KernelIterator kernel, KernelAccessor ka,
390 int kleft,
int kright,
391 int start = 0,
int stop = 0)
393 int w = std::distance( is, iend );
400 id += kright - start;
411 typedef typename PromoteTraits<
412 typename SrcAccessor::value_type,
413 typename KernelAccessor::value_type>::Promote SumType;
417 for(
int x=start; x<stop; ++x, ++is, ++id)
419 KernelIterator ik = kernel + kright;
420 SumType sum = NumericTraits<SumType>::zero();
422 SrcIterator iss = is + (-kright);
423 SrcIterator isend = is + (1 - kleft);
424 for(; iss != isend ; --ik, ++iss)
426 sum += ka(ik) * sa(iss);
429 da.set(detail::RequiresExplicitCast<
typename
430 DestAccessor::value_type>::cast(sum),
id);
576 doxygen_overloaded_function(template <...>
void convolveLine)
578 template <
class SrcIterator,
class SrcAccessor,
579 class DestIterator,
class DestAccessor,
580 class KernelIterator,
class KernelAccessor>
581 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
582 DestIterator
id, DestAccessor da,
583 KernelIterator ik, KernelAccessor ka,
584 int kleft,
int kright, BorderTreatmentMode border,
585 int start = 0,
int stop = 0)
587 typedef typename KernelAccessor::value_type KernelValue;
589 vigra_precondition(kleft <= 0,
590 "convolveLine(): kleft must be <= 0.\n");
591 vigra_precondition(kright >= 0,
592 "convolveLine(): kright must be >= 0.\n");
595 int w = std::distance( is, iend );
597 vigra_precondition(w >=
std::max(kright, -kleft) + 1,
598 "convolveLine(): kernel longer than line.\n");
601 vigra_precondition(0 <= start && start < stop && stop <= w,
602 "convolveLine(): invalid subrange (start, stop).\n");
606 case BORDER_TREATMENT_WRAP:
608 internalConvolveLineWrap(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
611 case BORDER_TREATMENT_AVOID:
613 internalConvolveLineAvoid(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
616 case BORDER_TREATMENT_REFLECT:
618 internalConvolveLineReflect(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
621 case BORDER_TREATMENT_REPEAT:
623 internalConvolveLineRepeat(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
626 case BORDER_TREATMENT_CLIP:
629 typedef typename KernelAccessor::value_type KT;
630 KT norm = NumericTraits<KT>::zero();
631 KernelIterator iik = ik + kleft;
632 for(
int i=kleft; i<=kright; ++i, ++iik)
635 vigra_precondition(norm != NumericTraits<KT>::zero(),
636 "convolveLine(): Norm of kernel must be != 0"
637 " in mode BORDER_TREATMENT_CLIP.\n");
639 internalConvolveLineClip(is, iend, sa,
id, da, ik, ka, kleft, kright, norm, start, stop);
644 vigra_precondition(0,
645 "convolveLine(): Unknown border treatment mode.\n");
650 template <
class SrcIterator,
class SrcAccessor,
651 class DestIterator,
class DestAccessor,
652 class KernelIterator,
class KernelAccessor>
654 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
655 pair<DestIterator, DestAccessor> dest,
656 tuple5<KernelIterator, KernelAccessor,
657 int,
int, BorderTreatmentMode> kernel,
658 int start = 0,
int stop = 0)
661 dest.first, dest.second,
662 kernel.first, kernel.second,
663 kernel.third, kernel.fourth, kernel.fifth, start, stop);
727 template <
class SrcIterator,
class SrcAccessor,
728 class DestIterator,
class DestAccessor,
729 class KernelIterator,
class KernelAccessor>
731 SrcIterator slowerright, SrcAccessor sa,
732 DestIterator dupperleft, DestAccessor da,
733 KernelIterator ik, KernelAccessor ka,
734 int kleft,
int kright, BorderTreatmentMode border)
736 typedef typename KernelAccessor::value_type KernelValue;
738 vigra_precondition(kleft <= 0,
739 "separableConvolveX(): kleft must be <= 0.\n");
740 vigra_precondition(kright >= 0,
741 "separableConvolveX(): kright must be >= 0.\n");
743 int w = slowerright.x - supperleft.x;
744 int h = slowerright.y - supperleft.y;
746 vigra_precondition(w >=
std::max(kright, -kleft) + 1,
747 "separableConvolveX(): kernel longer than line\n");
751 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
753 typename SrcIterator::row_iterator rs = supperleft.rowIterator();
754 typename DestIterator::row_iterator rd = dupperleft.rowIterator();
757 ik, ka, kleft, kright, border);
761 template <
class SrcIterator,
class SrcAccessor,
762 class DestIterator,
class DestAccessor,
763 class KernelIterator,
class KernelAccessor>
766 pair<DestIterator, DestAccessor> dest,
767 tuple5<KernelIterator, KernelAccessor,
768 int,
int, BorderTreatmentMode> kernel)
771 dest.first, dest.second,
772 kernel.first, kernel.second,
773 kernel.third, kernel.fourth, kernel.fifth);
839 template <
class SrcIterator,
class SrcAccessor,
840 class DestIterator,
class DestAccessor,
841 class KernelIterator,
class KernelAccessor>
843 SrcIterator slowerright, SrcAccessor sa,
844 DestIterator dupperleft, DestAccessor da,
845 KernelIterator ik, KernelAccessor ka,
846 int kleft,
int kright, BorderTreatmentMode border)
848 typedef typename KernelAccessor::value_type KernelValue;
850 vigra_precondition(kleft <= 0,
851 "separableConvolveY(): kleft must be <= 0.\n");
852 vigra_precondition(kright >= 0,
853 "separableConvolveY(): kright must be >= 0.\n");
855 int w = slowerright.x - supperleft.x;
856 int h = slowerright.y - supperleft.y;
858 vigra_precondition(h >=
std::max(kright, -kleft) + 1,
859 "separableConvolveY(): kernel longer than line\n");
863 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
865 typename SrcIterator::column_iterator cs = supperleft.columnIterator();
866 typename DestIterator::column_iterator cd = dupperleft.columnIterator();
869 ik, ka, kleft, kright, border);
873 template <
class SrcIterator,
class SrcAccessor,
874 class DestIterator,
class DestAccessor,
875 class KernelIterator,
class KernelAccessor>
878 pair<DestIterator, DestAccessor> dest,
879 tuple5<KernelIterator, KernelAccessor,
880 int,
int, BorderTreatmentMode> kernel)
883 dest.first, dest.second,
884 kernel.first, kernel.second,
885 kernel.third, kernel.fourth, kernel.fifth);
941 template <
class ARITHTYPE>
961 typedef typename InternalVector::iterator
Iterator;
965 typedef typename InternalVector::iterator
iterator;
982 : iter_(i), base_(i),
983 count_(count), sum_(count),
989 vigra_precondition(count_ == 1 || count_ == sum_,
990 "Kernel1D::initExplicitly(): "
991 "Wrong number of init values.");
1017 static value_type one() {
return NumericTraits<value_type>::one(); }
1027 border_treatment_(BORDER_TREATMENT_REFLECT),
1030 kernel_.push_back(norm_);
1036 : kernel_(k.kernel_),
1039 border_treatment_(k.border_treatment_),
1062 border_treatment_ = k.border_treatment_;
1064 kernel_ = k.kernel_;
1087 int size = right_ - left_ + 1;
1088 for(
unsigned int i=0; i<kernel_.
size(); ++i) kernel_[i] = v;
1089 norm_ = (double)size*v;
1091 return InitProxy(kernel_.
begin(),
size, norm_);
1307 this->
initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1335 this->
initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1363 this->
initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1398 vigra_precondition(a >= 0.0 && a <= 0.125,
1399 "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
1589 this->
initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
1629 vigra_precondition(left <= 0,
1630 "Kernel1D::initExplicitly(): left border must be <= 0.");
1631 vigra_precondition(right >= 0,
1632 "Kernel1D::initExplicitly(): right border must be >= 0.");
1637 kernel_.resize(right - left + 1);
1670 return kernel_[location -
left()];
1675 return kernel_[location -
left()];
1688 int size()
const {
return right_ - left_ + 1; }
1693 {
return border_treatment_; }
1698 { border_treatment_ = new_mode; }
1727 InternalVector kernel_;
1729 BorderTreatmentMode border_treatment_;
1733 template <
class ARITHTYPE>
1735 unsigned int derivativeOrder,
1738 typedef typename NumericTraits<value_type>::RealPromote TmpType;
1742 TmpType sum = NumericTraits<TmpType>::zero();
1744 if(derivativeOrder == 0)
1746 for(; k < kernel_.end(); ++k)
1753 unsigned int faculty = 1;
1754 for(
unsigned int i = 2; i <= derivativeOrder; ++i)
1756 for(
double x = left() + offset; k < kernel_.end(); ++x, ++k)
1758 sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x,
int(derivativeOrder)) / faculty);
1762 vigra_precondition(sum != NumericTraits<value_type>::zero(),
1763 "Kernel1D<ARITHTYPE>::normalize(): "
1764 "Cannot normalize a kernel with sum = 0");
1767 k = kernel_.begin();
1768 for(; k != kernel_.end(); ++k)
1778 template <
class ARITHTYPE>
1783 vigra_precondition(std_dev >= 0.0,
1784 "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
1785 vigra_precondition(windowRatio >= 0.0,
1786 "Kernel1D::initGaussian(): windowRatio must be >= 0.");
1794 if (windowRatio == 0.0)
1795 radius = (int)(3.0 * std_dev + 0.5);
1797 radius = (int)(windowRatio * std_dev + 0.5);
1802 kernel_.erase(kernel_.begin(), kernel_.end());
1803 kernel_.reserve(radius*2+1);
1805 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
1807 kernel_.push_back(gauss(x));
1814 kernel_.erase(kernel_.begin(), kernel_.end());
1815 kernel_.push_back(1.0);
1826 border_treatment_ = BORDER_TREATMENT_REFLECT;
1831 template <
class ARITHTYPE>
1835 vigra_precondition(std_dev >= 0.0,
1836 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
1841 int radius = (int)(3.0*std_dev + 0.5);
1845 double f = 2.0 / std_dev / std_dev;
1848 int maxIndex = (int)(2.0 * (radius + 5.0 *
VIGRA_CSTD::sqrt((
double)radius)) + 0.5);
1850 warray[maxIndex] = 0.0;
1851 warray[maxIndex-1] = 1.0;
1853 for(
int i = maxIndex-2; i >= radius; --i)
1855 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
1856 if(warray[i] > 1.0e40)
1858 warray[i+1] /= warray[i];
1866 warray[radius+1] = er * warray[radius+1] / warray[radius];
1867 warray[radius] = er;
1869 for(
int i = radius-1; i >= 0; --i)
1871 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
1875 double scale = norm / (2*er - warray[0]);
1877 initExplicitly(-radius, radius);
1880 for(
int i=0; i<=radius; ++i)
1882 c[i] = c[-i] = warray[i] * scale;
1887 kernel_.erase(kernel_.begin(), kernel_.end());
1888 kernel_.push_back(norm);
1896 border_treatment_ = BORDER_TREATMENT_REFLECT;
1901 template <
class ARITHTYPE>
1908 vigra_precondition(order >= 0,
1909 "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
1913 initGaussian(std_dev, norm, windowRatio);
1917 vigra_precondition(std_dev > 0.0,
1918 "Kernel1D::initGaussianDerivative(): "
1919 "Standard deviation must be > 0.");
1920 vigra_precondition(windowRatio >= 0.0,
1921 "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
1927 if(windowRatio == 0.0)
1928 radius = (int)(3.0 * std_dev + 0.5 * order + 0.5);
1930 radius = (int)(windowRatio * std_dev + 0.5);
1936 kernel_.reserve(radius*2+1);
1941 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
1943 kernel_.push_back(gauss(x));
1944 dc += kernel_[kernel_.size()-1];
1946 dc = ARITHTYPE(dc / (2.0*radius + 1.0));
1952 for(
unsigned int i=0; i < kernel_.size(); ++i)
1962 normalize(norm, order);
1968 border_treatment_ = BORDER_TREATMENT_REFLECT;
1973 template <
class ARITHTYPE>
1978 vigra_precondition(radius > 0,
1979 "Kernel1D::initBinomial(): Radius must be > 0.");
1983 typename InternalVector::iterator x = kernel_.begin() + radius;
1987 for(
int j=radius-1; j>=-radius; --j)
1989 x[j] = 0.5 * x[j+1];
1990 for(
int i=j+1; i<radius; ++i)
1992 x[i] = 0.5 * (x[i] + x[i+1]);
2002 border_treatment_ = BORDER_TREATMENT_REFLECT;
2007 template <
class ARITHTYPE>
2011 vigra_precondition(radius > 0,
2012 "Kernel1D::initAveraging(): Radius must be > 0.");
2015 double scale = 1.0 / (radius * 2 + 1);
2018 kernel_.erase(kernel_.begin(), kernel_.end());
2019 kernel_.reserve(radius*2+1);
2021 for(
int i=0; i<=radius*2+1; ++i)
2023 kernel_.push_back(scale * norm);
2031 border_treatment_ = BORDER_TREATMENT_CLIP;
2036 template <
class ARITHTYPE>
2040 kernel_.erase(kernel_.begin(), kernel_.end());
2043 kernel_.push_back(ARITHTYPE(0.5 * norm));
2044 kernel_.push_back(ARITHTYPE(0.0 * norm));
2045 kernel_.push_back(ARITHTYPE(-0.5 * norm));
2053 border_treatment_ = BORDER_TREATMENT_REFLECT;
2064 template <
class KernelIterator,
class KernelAccessor>
2066 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
2067 kernel1d(KernelIterator ik, KernelAccessor ka,
2068 int kleft,
int kright, BorderTreatmentMode border)
2071 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2072 ik, ka, kleft, kright, border);
2077 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2078 int, int, BorderTreatmentMode>
2079 kernel1d(Kernel1D<T>
const & k)
2083 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2084 int, int, BorderTreatmentMode>(
2087 k.left(), k.right(),
2088 k.borderTreatment());
2093 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2094 int, int, BorderTreatmentMode>
2095 kernel1d(Kernel1D<T>
const & k, BorderTreatmentMode border)
2099 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2100 int, int, BorderTreatmentMode>(
2103 k.left(), k.right(),
2110 #endif // VIGRA_SEPARABLECONVOLUTION_HXX