37 #ifndef VIGRA_SEPARABLECONVOLUTION_HXX
38 #define VIGRA_SEPARABLECONVOLUTION_HXX
42 #include "numerictraits.hxx"
43 #include "imageiteratoradapter.hxx"
44 #include "bordertreatment.hxx"
45 #include "gaussians.hxx"
46 #include "array_vector.hxx"
47 #include "multi_shape.hxx"
51 template <
class ARITHTYPE>
64 template <
class SrcIterator,
class SrcAccessor,
65 class DestIterator,
class DestAccessor,
66 class KernelIterator,
class KernelAccessor>
67 void internalConvolveLineOptimistic(SrcIterator is, SrcIterator iend, SrcAccessor sa,
68 DestIterator
id, DestAccessor da,
69 KernelIterator kernel, KernelAccessor ka,
70 int kleft,
int kright)
72 typedef typename PromoteTraits<
73 typename SrcAccessor::value_type,
74 typename KernelAccessor::value_type>::Promote SumType;
76 int w = std::distance( is, iend );
77 int kw = kright - kleft + 1;
78 for(
int x=0; x<w; ++x, ++is, ++id)
80 SrcIterator iss = is + (-kright);
81 KernelIterator ik = kernel + kright;
82 SumType
sum = NumericTraits<SumType>::zero();
84 for(
int k = 0; k < kw; ++k, --ik, ++iss)
86 sum += ka(ik) * sa(iss);
89 da.set(detail::RequiresExplicitCast<
typename
90 DestAccessor::value_type>::cast(sum),
id);
97 template <
class SrcIterator,
class SrcAccessor,
98 class DestIterator,
class DestAccessor>
100 copyLineWithBorderTreatment(SrcIterator is, SrcIterator iend, SrcAccessor sa,
101 DestIterator
id, DestAccessor da,
103 int kleft,
int kright,
104 BorderTreatmentMode borderTreatment)
106 int w = std::distance( is, iend );
107 int leftBorder = start - kright;
108 int rightBorder = stop - kleft;
109 int copyEnd = std::min(w, rightBorder);
113 switch(borderTreatment)
115 case BORDER_TREATMENT_WRAP:
117 for(; leftBorder<0; ++leftBorder, ++id)
118 da.set(sa(iend, leftBorder),
id);
121 case BORDER_TREATMENT_AVOID:
126 case BORDER_TREATMENT_REFLECT:
128 for(; leftBorder<0; ++leftBorder, ++id)
129 da.set(sa(is, -leftBorder),
id);
132 case BORDER_TREATMENT_REPEAT:
134 for(; leftBorder<0; ++leftBorder, ++id)
138 case BORDER_TREATMENT_CLIP:
140 vigra_precondition(
false,
141 "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
144 case BORDER_TREATMENT_ZEROPAD:
146 for(; leftBorder<0; ++leftBorder, ++id)
147 da.set(NumericTraits<typename DestAccessor::value_type>::zero(),
id);
152 vigra_precondition(
false,
153 "copyLineWithBorderTreatment(): Unknown border treatment mode.");
158 SrcIterator iss = is + leftBorder;
159 vigra_invariant( leftBorder < copyEnd,
160 "copyLineWithBorderTreatment(): assertion failed.");
161 for(; leftBorder<copyEnd; ++leftBorder, ++id, ++iss)
164 if(copyEnd < rightBorder)
166 switch(borderTreatment)
168 case BORDER_TREATMENT_WRAP:
170 for(; copyEnd<rightBorder; ++copyEnd, ++id, ++is)
174 case BORDER_TREATMENT_AVOID:
179 case BORDER_TREATMENT_REFLECT:
182 for(; copyEnd<rightBorder; ++copyEnd, ++id, --iss)
186 case BORDER_TREATMENT_REPEAT:
189 for(; copyEnd<rightBorder; ++copyEnd, ++id)
193 case BORDER_TREATMENT_CLIP:
195 vigra_precondition(
false,
196 "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
199 case BORDER_TREATMENT_ZEROPAD:
201 for(; copyEnd<rightBorder; ++copyEnd, ++id)
202 da.set(NumericTraits<typename DestAccessor::value_type>::zero(),
id);
207 vigra_precondition(
false,
208 "copyLineWithBorderTreatment(): Unknown border treatment mode.");
222 template <
class SrcIterator,
class SrcAccessor,
223 class DestIterator,
class DestAccessor,
224 class KernelIterator,
class KernelAccessor>
225 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
226 DestIterator
id, DestAccessor da,
227 KernelIterator kernel, KernelAccessor ka,
228 int kleft,
int kright,
229 int start = 0,
int stop = 0)
231 int w = std::distance( is, iend );
233 typedef typename PromoteTraits<
234 typename SrcAccessor::value_type,
235 typename KernelAccessor::value_type>::Promote SumType;
237 SrcIterator ibegin = is;
243 for(
int x=start; x<stop; ++x, ++is, ++id)
245 KernelIterator ik = kernel + kright;
246 SumType sum = NumericTraits<SumType>::zero();
251 SrcIterator iss = iend + x0;
253 for(; x0; ++x0, --ik, ++iss)
255 sum += ka(ik) * sa(iss);
261 SrcIterator isend = iend;
262 for(; iss != isend ; --ik, ++iss)
264 sum += ka(ik) * sa(iss);
267 int x0 = -kleft - w + x + 1;
270 for(; x0; --x0, --ik, ++iss)
272 sum += ka(ik) * sa(iss);
277 SrcIterator isend = is + (1 - kleft);
278 for(; iss != isend ; --ik, ++iss)
280 sum += ka(ik) * sa(iss);
284 else if(w-x <= -kleft)
286 SrcIterator iss = is + (-kright);
287 SrcIterator isend = iend;
288 for(; iss != isend ; --ik, ++iss)
290 sum += ka(ik) * sa(iss);
293 int x0 = -kleft - w + x + 1;
296 for(; x0; --x0, --ik, ++iss)
298 sum += ka(ik) * sa(iss);
303 SrcIterator iss = is - kright;
304 SrcIterator isend = is + (1 - kleft);
305 for(; iss != isend ; --ik, ++iss)
307 sum += ka(ik) * sa(iss);
311 da.set(detail::RequiresExplicitCast<
typename
312 DestAccessor::value_type>::cast(sum),
id);
322 template <
class SrcIterator,
class SrcAccessor,
323 class DestIterator,
class DestAccessor,
324 class KernelIterator,
class KernelAccessor,
326 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
327 DestIterator
id, DestAccessor da,
328 KernelIterator kernel, KernelAccessor ka,
329 int kleft,
int kright, Norm
norm,
330 int start = 0,
int stop = 0)
332 int w = std::distance( is, iend );
334 typedef typename PromoteTraits<
335 typename SrcAccessor::value_type,
336 typename KernelAccessor::value_type>::Promote SumType;
338 SrcIterator ibegin = is;
344 for(
int x=start; x<stop; ++x, ++is, ++id)
346 KernelIterator ik = kernel + kright;
347 SumType sum = NumericTraits<SumType>::zero();
352 Norm clipped = NumericTraits<Norm>::zero();
354 for(; x0; ++x0, --ik)
359 SrcIterator iss = ibegin;
362 SrcIterator isend = iend;
363 for(; iss != isend ; --ik, ++iss)
365 sum += ka(ik) * sa(iss);
368 int x0 = -kleft - w + x + 1;
370 for(; x0; --x0, --ik)
377 SrcIterator isend = is + (1 - kleft);
378 for(; iss != isend ; --ik, ++iss)
380 sum += ka(ik) * sa(iss);
384 sum = norm / (norm - clipped) * sum;
386 else if(w-x <= -kleft)
388 SrcIterator iss = is + (-kright);
389 SrcIterator isend = iend;
390 for(; iss != isend ; --ik, ++iss)
392 sum += ka(ik) * sa(iss);
395 Norm clipped = NumericTraits<Norm>::zero();
397 int x0 = -kleft - w + x + 1;
399 for(; x0; --x0, --ik)
404 sum = norm / (norm - clipped) * sum;
408 SrcIterator iss = is + (-kright);
409 SrcIterator isend = is + (1 - kleft);
410 for(; iss != isend ; --ik, ++iss)
412 sum += ka(ik) * sa(iss);
416 da.set(detail::RequiresExplicitCast<
typename
417 DestAccessor::value_type>::cast(sum),
id);
427 template <
class SrcIterator,
class SrcAccessor,
428 class DestIterator,
class DestAccessor,
429 class KernelIterator,
class KernelAccessor>
430 void internalConvolveLineZeropad(SrcIterator is, SrcIterator iend, SrcAccessor sa,
431 DestIterator
id, DestAccessor da,
432 KernelIterator kernel, KernelAccessor ka,
433 int kleft,
int kright,
434 int start = 0,
int stop = 0)
436 int w = std::distance( is, iend );
438 typedef typename PromoteTraits<
439 typename SrcAccessor::value_type,
440 typename KernelAccessor::value_type>::Promote SumType;
442 SrcIterator ibegin = is;
448 for(
int x=start; x<stop; ++x, ++is, ++id)
450 SumType sum = NumericTraits<SumType>::zero();
454 KernelIterator ik = kernel + x;
455 SrcIterator iss = ibegin;
459 SrcIterator isend = iend;
460 for(; iss != isend ; --ik, ++iss)
462 sum += ka(ik) * sa(iss);
467 SrcIterator isend = is + (1 - kleft);
468 for(; iss != isend ; --ik, ++iss)
470 sum += ka(ik) * sa(iss);
474 else if(w-x <= -kleft)
476 KernelIterator ik = kernel + kright;
477 SrcIterator iss = is + (-kright);
478 SrcIterator isend = iend;
479 for(; iss != isend ; --ik, ++iss)
481 sum += ka(ik) * sa(iss);
486 KernelIterator ik = kernel + kright;
487 SrcIterator iss = is + (-kright);
488 SrcIterator isend = is + (1 - kleft);
489 for(; iss != isend ; --ik, ++iss)
491 sum += ka(ik) * sa(iss);
495 da.set(detail::RequiresExplicitCast<
typename
496 DestAccessor::value_type>::cast(sum),
id);
506 template <
class SrcIterator,
class SrcAccessor,
507 class DestIterator,
class DestAccessor,
508 class KernelIterator,
class KernelAccessor>
509 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
510 DestIterator
id, DestAccessor da,
511 KernelIterator kernel, KernelAccessor ka,
512 int kleft,
int kright,
513 int start = 0,
int stop = 0)
515 int w = std::distance( is, iend );
517 typedef typename PromoteTraits<
518 typename SrcAccessor::value_type,
519 typename KernelAccessor::value_type>::Promote SumType;
521 SrcIterator ibegin = is;
527 for(
int x=start; x<stop; ++x, ++is, ++id)
529 KernelIterator ik = kernel + kright;
530 SumType sum = NumericTraits<SumType>::zero();
535 SrcIterator iss = ibegin - x0;
537 for(; x0; ++x0, --ik, --iss)
539 sum += ka(ik) * sa(iss);
544 SrcIterator isend = iend;
545 for(; iss != isend ; --ik, ++iss)
547 sum += ka(ik) * sa(iss);
550 int x0 = -kleft - w + x + 1;
553 for(; x0; --x0, --ik, --iss)
555 sum += ka(ik) * sa(iss);
560 SrcIterator isend = is + (1 - kleft);
561 for(; iss != isend ; --ik, ++iss)
563 sum += ka(ik) * sa(iss);
567 else if(w-x <= -kleft)
569 SrcIterator iss = is + (-kright);
570 SrcIterator isend = iend;
571 for(; iss != isend ; --ik, ++iss)
573 sum += ka(ik) * sa(iss);
576 int x0 = -kleft - w + x + 1;
579 for(; x0; --x0, --ik, --iss)
581 sum += ka(ik) * sa(iss);
586 SrcIterator iss = is + (-kright);
587 SrcIterator isend = is + (1 - kleft);
588 for(; iss != isend ; --ik, ++iss)
590 sum += ka(ik) * sa(iss);
594 da.set(detail::RequiresExplicitCast<
typename
595 DestAccessor::value_type>::cast(sum),
id);
605 template <
class SrcIterator,
class SrcAccessor,
606 class DestIterator,
class DestAccessor,
607 class KernelIterator,
class KernelAccessor>
608 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
609 DestIterator
id, DestAccessor da,
610 KernelIterator kernel, KernelAccessor ka,
611 int kleft,
int kright,
612 int start = 0,
int stop = 0)
614 int w = std::distance( is, iend );
616 typedef typename PromoteTraits<
617 typename SrcAccessor::value_type,
618 typename KernelAccessor::value_type>::Promote SumType;
620 SrcIterator ibegin = is;
626 for(
int x=start; x<stop; ++x, ++is, ++id)
628 KernelIterator ik = kernel + kright;
629 SumType sum = NumericTraits<SumType>::zero();
634 SrcIterator iss = ibegin;
636 for(; x0; ++x0, --ik)
638 sum += ka(ik) * sa(iss);
643 SrcIterator isend = iend;
644 for(; iss != isend ; --ik, ++iss)
646 sum += ka(ik) * sa(iss);
649 int x0 = -kleft - w + x + 1;
652 for(; x0; --x0, --ik)
654 sum += ka(ik) * sa(iss);
659 SrcIterator isend = is + (1 - kleft);
660 for(; iss != isend ; --ik, ++iss)
662 sum += ka(ik) * sa(iss);
666 else if(w-x <= -kleft)
668 SrcIterator iss = is + (-kright);
669 SrcIterator isend = iend;
670 for(; iss != isend ; --ik, ++iss)
672 sum += ka(ik) * sa(iss);
675 int x0 = -kleft - w + x + 1;
678 for(; x0; --x0, --ik)
680 sum += ka(ik) * sa(iss);
685 SrcIterator iss = is + (-kright);
686 SrcIterator isend = is + (1 - kleft);
687 for(; iss != isend ; --ik, ++iss)
689 sum += ka(ik) * sa(iss);
693 da.set(detail::RequiresExplicitCast<
typename
694 DestAccessor::value_type>::cast(sum),
id);
704 template <
class SrcIterator,
class SrcAccessor,
705 class DestIterator,
class DestAccessor,
706 class KernelIterator,
class KernelAccessor>
707 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
708 DestIterator
id, DestAccessor da,
709 KernelIterator kernel, KernelAccessor ka,
710 int kleft,
int kright,
711 int start = 0,
int stop = 0)
713 int w = std::distance( is, iend );
720 id += kright - start;
731 typedef typename PromoteTraits<
732 typename SrcAccessor::value_type,
733 typename KernelAccessor::value_type>::Promote SumType;
737 for(
int x=start; x<stop; ++x, ++is, ++id)
739 KernelIterator ik = kernel + kright;
740 SumType sum = NumericTraits<SumType>::zero();
742 SrcIterator iss = is + (-kright);
743 SrcIterator isend = is + (1 - kleft);
744 for(; iss != isend ; --ik, ++iss)
746 sum += ka(ik) * sa(iss);
749 da.set(detail::RequiresExplicitCast<
typename
750 DestAccessor::value_type>::cast(sum),
id);
896 template <
class SrcIterator,
class SrcAccessor,
897 class DestIterator,
class DestAccessor,
898 class KernelIterator,
class KernelAccessor>
899 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
900 DestIterator
id, DestAccessor da,
901 KernelIterator ik, KernelAccessor ka,
902 int kleft,
int kright, BorderTreatmentMode border,
903 int start = 0,
int stop = 0)
905 vigra_precondition(kleft <= 0,
906 "convolveLine(): kleft must be <= 0.\n");
907 vigra_precondition(kright >= 0,
908 "convolveLine(): kright must be >= 0.\n");
911 int w = std::distance( is, iend );
913 vigra_precondition(w >= std::max(kright, -kleft) + 1,
914 "convolveLine(): kernel longer than line.\n");
917 vigra_precondition(0 <= start && start < stop && stop <= w,
918 "convolveLine(): invalid subrange (start, stop).\n");
920 typedef typename PromoteTraits<
921 typename SrcAccessor::value_type,
922 typename KernelAccessor::value_type>::Promote SumType;
923 ArrayVector<SumType> a(iend - is);
926 case BORDER_TREATMENT_WRAP:
928 internalConvolveLineWrap(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
931 case BORDER_TREATMENT_AVOID:
933 internalConvolveLineAvoid(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
936 case BORDER_TREATMENT_REFLECT:
938 internalConvolveLineReflect(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
941 case BORDER_TREATMENT_REPEAT:
943 internalConvolveLineRepeat(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
946 case BORDER_TREATMENT_CLIP:
949 typedef typename KernelAccessor::value_type KT;
950 KT norm = NumericTraits<KT>::zero();
951 KernelIterator iik = ik + kleft;
952 for(
int i=kleft; i<=kright; ++i, ++iik)
955 vigra_precondition(norm != NumericTraits<KT>::zero(),
956 "convolveLine(): Norm of kernel must be != 0"
957 " in mode BORDER_TREATMENT_CLIP.\n");
959 internalConvolveLineClip(is, iend, sa,
id, da, ik, ka, kleft, kright, norm, start, stop);
962 case BORDER_TREATMENT_ZEROPAD:
964 internalConvolveLineZeropad(is, iend, sa,
id, da, ik, ka, kleft, kright, start, stop);
969 vigra_precondition(0,
970 "convolveLine(): Unknown border treatment mode.\n");
975 template <
class SrcIterator,
class SrcAccessor,
976 class DestIterator,
class DestAccessor,
977 class KernelIterator,
class KernelAccessor>
979 void convolveLine(triple<SrcIterator, SrcIterator, SrcAccessor> src,
980 pair<DestIterator, DestAccessor> dest,
981 tuple5<KernelIterator, KernelAccessor,
982 int,
int, BorderTreatmentMode> kernel,
983 int start = 0,
int stop = 0)
986 dest.first, dest.second,
987 kernel.first, kernel.second,
988 kernel.third, kernel.fourth, kernel.fifth, start, stop);
1085 template <
class SrcIterator,
class SrcAccessor,
1086 class DestIterator,
class DestAccessor,
1087 class KernelIterator,
class KernelAccessor>
1089 SrcIterator slowerright, SrcAccessor sa,
1090 DestIterator dupperleft, DestAccessor da,
1091 KernelIterator ik, KernelAccessor ka,
1092 int kleft,
int kright, BorderTreatmentMode border)
1094 vigra_precondition(kleft <= 0,
1095 "separableConvolveX(): kleft must be <= 0.\n");
1096 vigra_precondition(kright >= 0,
1097 "separableConvolveX(): kright must be >= 0.\n");
1099 int w = slowerright.x - supperleft.x;
1100 int h = slowerright.y - supperleft.y;
1102 vigra_precondition(w >= std::max(kright, -kleft) + 1,
1103 "separableConvolveX(): kernel longer than line\n");
1107 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1109 typename SrcIterator::row_iterator rs = supperleft.rowIterator();
1110 typename DestIterator::row_iterator rd = dupperleft.rowIterator();
1113 ik, ka, kleft, kright, border);
1117 template <
class SrcIterator,
class SrcAccessor,
1118 class DestIterator,
class DestAccessor,
1119 class KernelIterator,
class KernelAccessor>
1122 pair<DestIterator, DestAccessor> dest,
1123 tuple5<KernelIterator, KernelAccessor,
1124 int,
int, BorderTreatmentMode> kernel)
1127 dest.first, dest.second,
1128 kernel.first, kernel.second,
1129 kernel.third, kernel.fourth, kernel.fifth);
1132 template <
class T1,
class S1,
1137 MultiArrayView<2, T2, S2> dest,
1138 Kernel1D<T3>
const & kernel)
1140 vigra_precondition(src.shape() == dest.shape(),
1141 "separableConvolveX(): shape mismatch between input and output.");
1143 destImage(dest), kernel1d(kernel));
1240 template <
class SrcIterator,
class SrcAccessor,
1241 class DestIterator,
class DestAccessor,
1242 class KernelIterator,
class KernelAccessor>
1244 SrcIterator slowerright, SrcAccessor sa,
1245 DestIterator dupperleft, DestAccessor da,
1246 KernelIterator ik, KernelAccessor ka,
1247 int kleft,
int kright, BorderTreatmentMode border)
1249 vigra_precondition(kleft <= 0,
1250 "separableConvolveY(): kleft must be <= 0.\n");
1251 vigra_precondition(kright >= 0,
1252 "separableConvolveY(): kright must be >= 0.\n");
1254 int w = slowerright.x - supperleft.x;
1255 int h = slowerright.y - supperleft.y;
1257 vigra_precondition(h >= std::max(kright, -kleft) + 1,
1258 "separableConvolveY(): kernel longer than line\n");
1262 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1264 typename SrcIterator::column_iterator cs = supperleft.columnIterator();
1265 typename DestIterator::column_iterator cd = dupperleft.columnIterator();
1268 ik, ka, kleft, kright, border);
1272 template <
class SrcIterator,
class SrcAccessor,
1273 class DestIterator,
class DestAccessor,
1274 class KernelIterator,
class KernelAccessor>
1277 pair<DestIterator, DestAccessor> dest,
1278 tuple5<KernelIterator, KernelAccessor,
1279 int,
int, BorderTreatmentMode> kernel)
1282 dest.first, dest.second,
1283 kernel.first, kernel.second,
1284 kernel.third, kernel.fourth, kernel.fifth);
1287 template <
class T1,
class S1,
1292 MultiArrayView<2, T2, S2> dest,
1293 Kernel1D<T3>
const & kernel)
1295 vigra_precondition(src.shape() == dest.shape(),
1296 "separableConvolveY(): shape mismatch between input and output.");
1298 destImage(dest), kernel1d(kernel));
1365 template <
class ARITHTYPE =
double>
1369 typedef ArrayVector<ARITHTYPE> InternalVector;
1406 : iter_(i), base_(i),
1407 count_(count), sum_(count),
1413 throw(PreconditionViolation)
1414 #elif _MSC_VER >= 1900
1418 vigra_precondition(count_ == 1 || count_ == sum_,
1419 "Kernel1D::initExplicitly(): "
1420 "Wrong number of init values.");
1445 static value_type one() {
return NumericTraits<value_type>::one(); }
1455 border_treatment_(BORDER_TREATMENT_REFLECT),
1458 kernel_.push_back(norm_);
1464 : kernel_(k.kernel_),
1467 border_treatment_(k.border_treatment_),
1478 border_treatment_(k.borderTreatment()),
1490 border_treatment_ = k.border_treatment_;
1492 kernel_ = k.kernel_;
1515 int size = right_ - left_ + 1;
1516 for(
unsigned int i=0; i<kernel_.
size(); ++i) kernel_[i] = v;
1517 norm_ = (double)size*v;
1519 return InitProxy(kernel_.
begin(),
size, norm_);
1735 this->
initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1763 this->
initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1791 this->
initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1826 vigra_precondition(a >= 0.0 && a <= 0.125,
1827 "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
2062 this->
initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
2102 vigra_precondition(left <= 0,
2103 "Kernel1D::initExplicitly(): left border must be <= 0.");
2104 vigra_precondition(right >= 0,
2105 "Kernel1D::initExplicitly(): right border must be >= 0.");
2110 kernel_.resize(right - left + 1);
2143 return kernel_[location -
left()];
2148 return kernel_[location -
left()];
2161 int size()
const {
return right_ - left_ + 1; }
2166 {
return border_treatment_; }
2171 { border_treatment_ = new_mode; }
2200 InternalVector kernel_;
2202 BorderTreatmentMode border_treatment_;
2206 template <
class ARITHTYPE>
2208 unsigned int derivativeOrder,
2211 typedef typename NumericTraits<value_type>::RealPromote TmpType;
2215 TmpType sum = NumericTraits<TmpType>::zero();
2217 if(derivativeOrder == 0)
2219 for(; k < kernel_.end(); ++k)
2226 unsigned int faculty = 1;
2227 for(
unsigned int i = 2; i <= derivativeOrder; ++i)
2229 for(
double x = left() + offset; k < kernel_.end(); ++x, ++k)
2231 sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x,
int(derivativeOrder)) / faculty);
2235 vigra_precondition(sum != NumericTraits<value_type>::zero(),
2236 "Kernel1D<ARITHTYPE>::normalize(): "
2237 "Cannot normalize a kernel with sum = 0");
2240 k = kernel_.begin();
2241 for(; k != kernel_.end(); ++k)
2251 template <
class ARITHTYPE>
2257 vigra_precondition(std_dev >= 0.0,
2258 "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
2259 vigra_precondition(windowRatio >= 0.0,
2260 "Kernel1D::initGaussian(): windowRatio must be >= 0.");
2268 if (windowRatio == 0.0)
2269 radius = (int)(3.0 * std_dev + 0.5);
2271 radius = (int)(windowRatio * std_dev + 0.5);
2276 kernel_.erase(kernel_.begin(), kernel_.end());
2277 kernel_.reserve(radius*2+1);
2279 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2281 kernel_.push_back(gauss(x));
2288 kernel_.erase(kernel_.begin(), kernel_.end());
2289 kernel_.push_back(1.0);
2300 border_treatment_ = BORDER_TREATMENT_REFLECT;
2305 template <
class ARITHTYPE>
2310 vigra_precondition(std_dev >= 0.0,
2311 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
2316 int radius = (int)(3.0*std_dev + 0.5);
2320 double f = 2.0 / std_dev / std_dev;
2323 int maxIndex = (int)(2.0 * (radius + 5.0 *
VIGRA_CSTD::sqrt((
double)radius)) + 0.5);
2325 warray[maxIndex] = 0.0;
2326 warray[maxIndex-1] = 1.0;
2328 for(
int i = maxIndex-2; i >= radius; --i)
2330 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2331 if(warray[i] > 1.0e40)
2333 warray[i+1] /= warray[i];
2341 warray[radius+1] = er * warray[radius+1] / warray[radius];
2342 warray[radius] = er;
2344 for(
int i = radius-1; i >= 0; --i)
2346 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2350 double scale = norm / (2*er - warray[0]);
2352 initExplicitly(-radius, radius);
2355 for(
int i=0; i<=radius; ++i)
2357 c[i] = c[-i] = warray[i] * scale;
2362 kernel_.erase(kernel_.begin(), kernel_.end());
2363 kernel_.push_back(norm);
2371 border_treatment_ = BORDER_TREATMENT_REFLECT;
2376 template <
class ARITHTYPE>
2383 vigra_precondition(order >= 0,
2384 "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
2388 initGaussian(std_dev, norm, windowRatio);
2392 vigra_precondition(std_dev > 0.0,
2393 "Kernel1D::initGaussianDerivative(): "
2394 "Standard deviation must be > 0.");
2395 vigra_precondition(windowRatio >= 0.0,
2396 "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
2402 if(windowRatio == 0.0)
2403 radius = (int)((3.0 + 0.5 * order) * std_dev + 0.5);
2405 radius = (int)(windowRatio * std_dev + 0.5);
2411 kernel_.reserve(radius*2+1);
2416 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2418 kernel_.push_back(gauss(x));
2419 dc += kernel_[kernel_.size()-1];
2421 dc = ARITHTYPE(dc / (2.0*radius + 1.0));
2427 for(
unsigned int i=0; i < kernel_.size(); ++i)
2437 normalize(norm, order);
2443 border_treatment_ = BORDER_TREATMENT_REFLECT;
2448 template <
class ARITHTYPE>
2453 vigra_precondition(radius > 0,
2454 "Kernel1D::initBinomial(): Radius must be > 0.");
2458 typename InternalVector::iterator x = kernel_.begin() + radius;
2462 for(
int j=radius-1; j>=-radius; --j)
2464 x[j] = 0.5 * x[j+1];
2465 for(
int i=j+1; i<radius; ++i)
2467 x[i] = 0.5 * (x[i] + x[i+1]);
2477 border_treatment_ = BORDER_TREATMENT_REFLECT;
2482 template <
class ARITHTYPE>
2487 vigra_precondition(radius > 0,
2488 "Kernel1D::initAveraging(): Radius must be > 0.");
2491 double scale = 1.0 / (radius * 2 + 1);
2494 kernel_.erase(kernel_.begin(), kernel_.end());
2495 kernel_.reserve(radius*2+1);
2497 for(
int i=0; i<=radius*2+1; ++i)
2499 kernel_.push_back(scale * norm);
2507 border_treatment_ = BORDER_TREATMENT_CLIP;
2512 template <
class ARITHTYPE>
2516 kernel_.erase(kernel_.begin(), kernel_.end());
2519 kernel_.push_back(ARITHTYPE(0.5 * norm));
2520 kernel_.push_back(ARITHTYPE(0.0 * norm));
2521 kernel_.push_back(ARITHTYPE(-0.5 * norm));
2529 border_treatment_ = BORDER_TREATMENT_REFLECT;
2540 template <
class KernelIterator,
class KernelAccessor>
2542 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
2543 kernel1d(KernelIterator ik, KernelAccessor ka,
2544 int kleft,
int kright, BorderTreatmentMode border)
2547 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2548 ik, ka, kleft, kright, border);
2553 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2554 int, int, BorderTreatmentMode>
2555 kernel1d(Kernel1D<T>
const & k)
2559 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2560 int, int, BorderTreatmentMode>(
2563 k.left(), k.right(),
2564 k.borderTreatment());
2569 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2570 int, int, BorderTreatmentMode>
2571 kernel1d(Kernel1D<T>
const & k, BorderTreatmentMode border)
2575 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2576 int, int, BorderTreatmentMode>(
2579 k.left(), k.right(),
2586 #endif // VIGRA_SEPARABLECONVOLUTION_HXX
StandardAccessor< ARITHTYPE > Accessor
Definition: separableconvolution.hxx:1397
InternalVector::const_reference const_reference
Definition: separableconvolution.hxx:1381
void normalize()
Definition: separableconvolution.hxx:2186
StandardConstAccessor< ARITHTYPE > ConstAccessor
Definition: separableconvolution.hxx:1401
value_type norm() const
Definition: separableconvolution.hxx:2175
void initAveraging(int radius, value_type norm)
Definition: separableconvolution.hxx:2484
void initOptimalFirstDerivative5()
Definition: separableconvolution.hxx:2031
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:52
void initBinomial(int radius, value_type norm)
Definition: separableconvolution.hxx:2450
Kernel1D(Kernel1D const &k)
Definition: separableconvolution.hxx:1463
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2253
InternalVector::value_type value_type
Definition: separableconvolution.hxx:1373
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void initGaussianDerivative(double std_dev, int order)
Definition: separableconvolution.hxx:1626
const_iterator begin() const
Definition: array_vector.hxx:223
InternalVector::iterator iterator
Definition: separableconvolution.hxx:1389
void initGaussian(double std_dev)
Definition: separableconvolution.hxx:1556
void initOptimalSecondDerivative5()
Definition: separableconvolution.hxx:2060
void initDiscreteGaussian(double std_dev)
Definition: separableconvolution.hxx:1584
void initSymmetricDifference()
Definition: separableconvolution.hxx:1978
Kernel1D()
Definition: separableconvolution.hxx:1451
void convolveLine(...)
Performs a 1-dimensional convolution of the source signal using the given kernel. ...
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
void setBorderTreatment(BorderTreatmentMode new_mode)
Definition: separableconvolution.hxx:2170
void initForwardDifference()
Definition: separableconvolution.hxx:1931
int left() const
Definition: separableconvolution.hxx:2153
Kernel1D & operator=(Kernel1D const &k)
Definition: separableconvolution.hxx:1484
int size() const
Definition: separableconvolution.hxx:2161
void initSecondDifference3()
Definition: separableconvolution.hxx:1999
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
Kernel1D & initExplicitly(int left, int right)
Definition: separableconvolution.hxx:2100
Definition: gaussians.hxx:63
void initDiscreteGaussian(double std_dev, value_type norm)
Definition: separableconvolution.hxx:2307
void initOptimalFirstDerivativeSmoothing3()
Definition: separableconvolution.hxx:1679
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
void initBurtFilter(double a=0.04785)
Definition: separableconvolution.hxx:1824
void initSymmetricGradient()
Definition: separableconvolution.hxx:1908
Accessor accessor()
Definition: separableconvolution.hxx:2197
void initSymmetricGradient(value_type norm)
Definition: separableconvolution.hxx:1898
Kernel1D(Kernel1D< U > const &k)
Definition: separableconvolution.hxx:1474
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void initOptimalSecondDerivativeSmoothing5()
Definition: separableconvolution.hxx:1789
InternalVector::reference reference
Definition: separableconvolution.hxx:1377
void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2378
void initOptimalSmoothing5()
Definition: separableconvolution.hxx:1733
Encapsulate read access to the values an iterator points to.
Definition: accessor.hxx:269
~Kernel1D()
Definition: separableconvolution.hxx:1524
BorderTreatmentMode borderTreatment() const
Definition: separableconvolution.hxx:2165
InternalVector::iterator Iterator
Definition: separableconvolution.hxx:1385
int right() const
Definition: separableconvolution.hxx:2157
InternalVector::const_iterator const_iterator
Definition: separableconvolution.hxx:1393
iterator center()
Definition: separableconvolution.hxx:2123
void initOptimalFirstDerivativeSmoothing5()
Definition: separableconvolution.hxx:1761
size_type size() const
Definition: array_vector.hxx:358
void initOptimalSecondDerivativeSmoothing3()
Definition: separableconvolution.hxx:1707
void initAveraging(int radius)
Definition: separableconvolution.hxx:1879
void initBinomial(int radius)
Definition: separableconvolution.hxx:1853
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:133
InitProxy operator=(value_type const &v)
Definition: separableconvolution.hxx:1513
void separableConvolveY(...)
Performs a 1 dimensional convolution in y direction.
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
reference operator[](int location)
Definition: separableconvolution.hxx:2141
ConstAccessor accessor() const
Definition: separableconvolution.hxx:2193
void initBackwardDifference()
Definition: separableconvolution.hxx:1955
void initOptimalSmoothing3()
Definition: separableconvolution.hxx:1651