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"
60 template <
class SrcIterator,
class SrcAccessor,
61 class DestIterator,
class DestAccessor,
62 class KernelIterator,
class KernelAccessor>
63 void internalConvolveLineOptimistic(SrcIterator is, SrcIterator iend, SrcAccessor sa,
64 DestIterator
id, DestAccessor da,
65 KernelIterator kernel, KernelAccessor ka,
66 int kleft,
int kright)
68 typedef typename PromoteTraits<
69 typename SrcAccessor::value_type,
70 typename KernelAccessor::value_type>::Promote SumType;
72 int w = std::distance( is, iend );
73 int kw = kright - kleft + 1;
74 for(
int x=0; x<w; ++x, ++is, ++id)
76 SrcIterator iss = is + (-kright);
77 KernelIterator ik = kernel + kright;
78 SumType
sum = NumericTraits<SumType>::zero();
80 for(
int k = 0; k < kw; ++k, --ik, ++iss)
82 sum += ka(ik) * sa(iss);
85 da.set(detail::RequiresExplicitCast<
typename
86 DestAccessor::value_type>::cast(sum),
id);
93 template <
class SrcIterator,
class SrcAccessor,
94 class DestIterator,
class DestAccessor>
96 copyLineWithBorderTreatment(SrcIterator is, SrcIterator iend, SrcAccessor sa,
97 DestIterator
id, DestAccessor da,
99 int kleft,
int kright,
100 BorderTreatmentMode borderTreatment)
102 int w = std::distance( is, iend );
103 int leftBorder = start - kright;
104 int rightBorder = stop - kleft;
105 int copyEnd = std::min(w, rightBorder);
109 switch(borderTreatment)
111 case BORDER_TREATMENT_WRAP:
113 for(; leftBorder<0; ++leftBorder, ++id)
114 da.set(sa(iend, leftBorder),
id);
117 case BORDER_TREATMENT_AVOID:
122 case BORDER_TREATMENT_REFLECT:
124 for(; leftBorder<0; ++leftBorder, ++id)
125 da.set(sa(is, -leftBorder),
id);
128 case BORDER_TREATMENT_REPEAT:
130 for(; leftBorder<0; ++leftBorder, ++id)
134 case BORDER_TREATMENT_CLIP:
136 vigra_precondition(
false,
137 "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
140 case BORDER_TREATMENT_ZEROPAD:
142 for(; leftBorder<0; ++leftBorder, ++id)
143 da.set(NumericTraits<typename DestAccessor::value_type>::zero(),
id);
148 vigra_precondition(
false,
149 "copyLineWithBorderTreatment(): Unknown border treatment mode.");
154 SrcIterator iss = is + leftBorder;
155 vigra_invariant( leftBorder < copyEnd,
156 "copyLineWithBorderTreatment(): assertion failed.");
157 for(; leftBorder<copyEnd; ++leftBorder, ++id, ++iss)
160 if(copyEnd < rightBorder)
162 switch(borderTreatment)
164 case BORDER_TREATMENT_WRAP:
166 for(; copyEnd<rightBorder; ++copyEnd, ++id, ++is)
170 case BORDER_TREATMENT_AVOID:
175 case BORDER_TREATMENT_REFLECT:
178 for(; copyEnd<rightBorder; ++copyEnd, ++id, --iss)
182 case BORDER_TREATMENT_REPEAT:
185 for(; copyEnd<rightBorder; ++copyEnd, ++id)
189 case BORDER_TREATMENT_CLIP:
191 vigra_precondition(
false,
192 "copyLineWithBorderTreatment() internal error: not applicable to BORDER_TREATMENT_CLIP.");
195 case BORDER_TREATMENT_ZEROPAD:
197 for(; copyEnd<rightBorder; ++copyEnd, ++id)
198 da.set(NumericTraits<typename DestAccessor::value_type>::zero(),
id);
203 vigra_precondition(
false,
204 "copyLineWithBorderTreatment(): Unknown border treatment mode.");
218 template <
class SrcIterator,
class SrcAccessor,
219 class DestIterator,
class DestAccessor,
220 class KernelIterator,
class KernelAccessor>
221 void internalConvolveLineWrap(SrcIterator is, SrcIterator iend, SrcAccessor sa,
222 DestIterator
id, DestAccessor da,
223 KernelIterator kernel, KernelAccessor ka,
224 int kleft,
int kright,
225 int start = 0,
int stop = 0)
227 int w = std::distance( is, iend );
229 typedef typename PromoteTraits<
230 typename SrcAccessor::value_type,
231 typename KernelAccessor::value_type>::Promote SumType;
233 SrcIterator ibegin = is;
239 for(
int x=start; x<stop; ++x, ++is, ++id)
241 KernelIterator ik = kernel + kright;
242 SumType sum = NumericTraits<SumType>::zero();
247 SrcIterator iss = iend + x0;
249 for(; x0; ++x0, --ik, ++iss)
251 sum += ka(ik) * sa(iss);
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 isend = is + (1 - kleft);
274 for(; iss != isend ; --ik, ++iss)
276 sum += ka(ik) * sa(iss);
280 else if(w-x <= -kleft)
282 SrcIterator iss = is + (-kright);
283 SrcIterator isend = iend;
284 for(; iss != isend ; --ik, ++iss)
286 sum += ka(ik) * sa(iss);
289 int x0 = -kleft - w + x + 1;
292 for(; x0; --x0, --ik, ++iss)
294 sum += ka(ik) * sa(iss);
299 SrcIterator iss = is - kright;
300 SrcIterator isend = is + (1 - kleft);
301 for(; iss != isend ; --ik, ++iss)
303 sum += ka(ik) * sa(iss);
307 da.set(detail::RequiresExplicitCast<
typename
308 DestAccessor::value_type>::cast(sum),
id);
318 template <
class SrcIterator,
class SrcAccessor,
319 class DestIterator,
class DestAccessor,
320 class KernelIterator,
class KernelAccessor,
322 void internalConvolveLineClip(SrcIterator is, SrcIterator iend, SrcAccessor sa,
323 DestIterator
id, DestAccessor da,
324 KernelIterator kernel, KernelAccessor ka,
325 int kleft,
int kright, Norm
norm,
326 int start = 0,
int stop = 0)
328 int w = std::distance( is, iend );
330 typedef typename PromoteTraits<
331 typename SrcAccessor::value_type,
332 typename KernelAccessor::value_type>::Promote SumType;
334 SrcIterator ibegin = is;
340 for(
int x=start; x<stop; ++x, ++is, ++id)
342 KernelIterator ik = kernel + kright;
343 SumType sum = NumericTraits<SumType>::zero();
348 Norm clipped = NumericTraits<Norm>::zero();
350 for(; x0; ++x0, --ik)
355 SrcIterator iss = ibegin;
358 SrcIterator isend = iend;
359 for(; iss != isend ; --ik, ++iss)
361 sum += ka(ik) * sa(iss);
364 int x0 = -kleft - w + x + 1;
366 for(; x0; --x0, --ik)
373 SrcIterator isend = is + (1 - kleft);
374 for(; iss != isend ; --ik, ++iss)
376 sum += ka(ik) * sa(iss);
380 sum = norm / (norm - clipped) * sum;
382 else if(w-x <= -kleft)
384 SrcIterator iss = is + (-kright);
385 SrcIterator isend = iend;
386 for(; iss != isend ; --ik, ++iss)
388 sum += ka(ik) * sa(iss);
391 Norm clipped = NumericTraits<Norm>::zero();
393 int x0 = -kleft - w + x + 1;
395 for(; x0; --x0, --ik)
400 sum = norm / (norm - clipped) * sum;
404 SrcIterator iss = is + (-kright);
405 SrcIterator isend = is + (1 - kleft);
406 for(; iss != isend ; --ik, ++iss)
408 sum += ka(ik) * sa(iss);
412 da.set(detail::RequiresExplicitCast<
typename
413 DestAccessor::value_type>::cast(sum),
id);
423 template <
class SrcIterator,
class SrcAccessor,
424 class DestIterator,
class DestAccessor,
425 class KernelIterator,
class KernelAccessor>
426 void internalConvolveLineZeropad(SrcIterator is, SrcIterator iend, SrcAccessor sa,
427 DestIterator
id, DestAccessor da,
428 KernelIterator kernel, KernelAccessor ka,
429 int kleft,
int kright,
430 int start = 0,
int stop = 0)
432 int w = std::distance( is, iend );
434 typedef typename PromoteTraits<
435 typename SrcAccessor::value_type,
436 typename KernelAccessor::value_type>::Promote SumType;
438 SrcIterator ibegin = is;
444 for(
int x=start; x<stop; ++x, ++is, ++id)
446 SumType sum = NumericTraits<SumType>::zero();
450 KernelIterator ik = kernel + x;
451 SrcIterator iss = ibegin;
455 SrcIterator isend = iend;
456 for(; iss != isend ; --ik, ++iss)
458 sum += ka(ik) * sa(iss);
463 SrcIterator isend = is + (1 - kleft);
464 for(; iss != isend ; --ik, ++iss)
466 sum += ka(ik) * sa(iss);
470 else if(w-x <= -kleft)
472 KernelIterator ik = kernel + kright;
473 SrcIterator iss = is + (-kright);
474 SrcIterator isend = iend;
475 for(; iss != isend ; --ik, ++iss)
477 sum += ka(ik) * sa(iss);
482 KernelIterator ik = kernel + kright;
483 SrcIterator iss = is + (-kright);
484 SrcIterator isend = is + (1 - kleft);
485 for(; iss != isend ; --ik, ++iss)
487 sum += ka(ik) * sa(iss);
491 da.set(detail::RequiresExplicitCast<
typename
492 DestAccessor::value_type>::cast(sum),
id);
502 template <
class SrcIterator,
class SrcAccessor,
503 class DestIterator,
class DestAccessor,
504 class KernelIterator,
class KernelAccessor>
505 void internalConvolveLineReflect(SrcIterator is, SrcIterator iend, SrcAccessor sa,
506 DestIterator
id, DestAccessor da,
507 KernelIterator kernel, KernelAccessor ka,
508 int kleft,
int kright,
509 int start = 0,
int stop = 0)
511 int w = std::distance( is, iend );
513 typedef typename PromoteTraits<
514 typename SrcAccessor::value_type,
515 typename KernelAccessor::value_type>::Promote SumType;
517 SrcIterator ibegin = is;
523 for(
int x=start; x<stop; ++x, ++is, ++id)
525 KernelIterator ik = kernel + kright;
526 SumType sum = NumericTraits<SumType>::zero();
531 SrcIterator iss = ibegin - x0;
533 for(; x0; ++x0, --ik, --iss)
535 sum += ka(ik) * sa(iss);
540 SrcIterator isend = iend;
541 for(; iss != isend ; --ik, ++iss)
543 sum += ka(ik) * sa(iss);
546 int x0 = -kleft - w + x + 1;
549 for(; x0; --x0, --ik, --iss)
551 sum += ka(ik) * sa(iss);
556 SrcIterator isend = is + (1 - kleft);
557 for(; iss != isend ; --ik, ++iss)
559 sum += ka(ik) * sa(iss);
563 else if(w-x <= -kleft)
565 SrcIterator iss = is + (-kright);
566 SrcIterator isend = iend;
567 for(; iss != isend ; --ik, ++iss)
569 sum += ka(ik) * sa(iss);
572 int x0 = -kleft - w + x + 1;
575 for(; x0; --x0, --ik, --iss)
577 sum += ka(ik) * sa(iss);
582 SrcIterator iss = is + (-kright);
583 SrcIterator isend = is + (1 - kleft);
584 for(; iss != isend ; --ik, ++iss)
586 sum += ka(ik) * sa(iss);
590 da.set(detail::RequiresExplicitCast<
typename
591 DestAccessor::value_type>::cast(sum),
id);
601 template <
class SrcIterator,
class SrcAccessor,
602 class DestIterator,
class DestAccessor,
603 class KernelIterator,
class KernelAccessor>
604 void internalConvolveLineRepeat(SrcIterator is, SrcIterator iend, SrcAccessor sa,
605 DestIterator
id, DestAccessor da,
606 KernelIterator kernel, KernelAccessor ka,
607 int kleft,
int kright,
608 int start = 0,
int stop = 0)
610 int w = std::distance( is, iend );
612 typedef typename PromoteTraits<
613 typename SrcAccessor::value_type,
614 typename KernelAccessor::value_type>::Promote SumType;
616 SrcIterator ibegin = is;
622 for(
int x=start; x<stop; ++x, ++is, ++id)
624 KernelIterator ik = kernel + kright;
625 SumType sum = NumericTraits<SumType>::zero();
630 SrcIterator iss = ibegin;
632 for(; x0; ++x0, --ik)
634 sum += ka(ik) * sa(iss);
639 SrcIterator isend = iend;
640 for(; iss != isend ; --ik, ++iss)
642 sum += ka(ik) * sa(iss);
645 int x0 = -kleft - w + x + 1;
648 for(; x0; --x0, --ik)
650 sum += ka(ik) * sa(iss);
655 SrcIterator isend = is + (1 - kleft);
656 for(; iss != isend ; --ik, ++iss)
658 sum += ka(ik) * sa(iss);
662 else if(w-x <= -kleft)
664 SrcIterator iss = is + (-kright);
665 SrcIterator isend = iend;
666 for(; iss != isend ; --ik, ++iss)
668 sum += ka(ik) * sa(iss);
671 int x0 = -kleft - w + x + 1;
674 for(; x0; --x0, --ik)
676 sum += ka(ik) * sa(iss);
681 SrcIterator iss = is + (-kright);
682 SrcIterator isend = is + (1 - kleft);
683 for(; iss != isend ; --ik, ++iss)
685 sum += ka(ik) * sa(iss);
689 da.set(detail::RequiresExplicitCast<
typename
690 DestAccessor::value_type>::cast(sum),
id);
700 template <
class SrcIterator,
class SrcAccessor,
701 class DestIterator,
class DestAccessor,
702 class KernelIterator,
class KernelAccessor>
703 void internalConvolveLineAvoid(SrcIterator is, SrcIterator iend, SrcAccessor sa,
704 DestIterator
id, DestAccessor da,
705 KernelIterator kernel, KernelAccessor ka,
706 int kleft,
int kright,
707 int start = 0,
int stop = 0)
709 int w = std::distance( is, iend );
716 id += kright - start;
727 typedef typename PromoteTraits<
728 typename SrcAccessor::value_type,
729 typename KernelAccessor::value_type>::Promote SumType;
733 for(
int x=start; x<stop; ++x, ++is, ++id)
735 KernelIterator ik = kernel + kright;
736 SumType sum = NumericTraits<SumType>::zero();
738 SrcIterator iss = is + (-kright);
739 SrcIterator isend = is + (1 - kleft);
740 for(; iss != isend ; --ik, ++iss)
742 sum += ka(ik) * sa(iss);
745 da.set(detail::RequiresExplicitCast<
typename
746 DestAccessor::value_type>::cast(sum),
id);
892 doxygen_overloaded_function(template <...>
void convolveLine)
894 template <
class SrcIterator,
class SrcAccessor,
895 class DestIterator,
class DestAccessor,
896 class KernelIterator,
class KernelAccessor>
897 void convolveLine(SrcIterator is, SrcIterator iend, SrcAccessor sa,
898 DestIterator
id, DestAccessor da,
899 KernelIterator ik, KernelAccessor ka,
900 int kleft,
int kright, BorderTreatmentMode border,
901 int start = 0,
int stop = 0)
903 typedef typename KernelAccessor::value_type KernelValue;
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);
1052 template <
class SrcIterator,
class SrcAccessor,
1053 class DestIterator,
class DestAccessor,
1054 class KernelIterator,
class KernelAccessor>
1056 SrcIterator slowerright, SrcAccessor sa,
1057 DestIterator dupperleft, DestAccessor da,
1058 KernelIterator ik, KernelAccessor ka,
1059 int kleft,
int kright, BorderTreatmentMode border)
1061 typedef typename KernelAccessor::value_type KernelValue;
1063 vigra_precondition(kleft <= 0,
1064 "separableConvolveX(): kleft must be <= 0.\n");
1065 vigra_precondition(kright >= 0,
1066 "separableConvolveX(): kright must be >= 0.\n");
1068 int w = slowerright.x - supperleft.x;
1069 int h = slowerright.y - supperleft.y;
1071 vigra_precondition(w >= std::max(kright, -kleft) + 1,
1072 "separableConvolveX(): kernel longer than line\n");
1076 for(y=0; y<h; ++y, ++supperleft.y, ++dupperleft.y)
1078 typename SrcIterator::row_iterator rs = supperleft.rowIterator();
1079 typename DestIterator::row_iterator rd = dupperleft.rowIterator();
1082 ik, ka, kleft, kright, border);
1086 template <
class SrcIterator,
class SrcAccessor,
1087 class DestIterator,
class DestAccessor,
1088 class KernelIterator,
class KernelAccessor>
1091 pair<DestIterator, DestAccessor> dest,
1092 tuple5<KernelIterator, KernelAccessor,
1093 int,
int, BorderTreatmentMode> kernel)
1096 dest.first, dest.second,
1097 kernel.first, kernel.second,
1098 kernel.third, kernel.fourth, kernel.fifth);
1164 template <
class SrcIterator,
class SrcAccessor,
1165 class DestIterator,
class DestAccessor,
1166 class KernelIterator,
class KernelAccessor>
1168 SrcIterator slowerright, SrcAccessor sa,
1169 DestIterator dupperleft, DestAccessor da,
1170 KernelIterator ik, KernelAccessor ka,
1171 int kleft,
int kright, BorderTreatmentMode border)
1173 typedef typename KernelAccessor::value_type KernelValue;
1175 vigra_precondition(kleft <= 0,
1176 "separableConvolveY(): kleft must be <= 0.\n");
1177 vigra_precondition(kright >= 0,
1178 "separableConvolveY(): kright must be >= 0.\n");
1180 int w = slowerright.x - supperleft.x;
1181 int h = slowerright.y - supperleft.y;
1183 vigra_precondition(h >= std::max(kright, -kleft) + 1,
1184 "separableConvolveY(): kernel longer than line\n");
1188 for(x=0; x<w; ++x, ++supperleft.x, ++dupperleft.x)
1190 typename SrcIterator::column_iterator cs = supperleft.columnIterator();
1191 typename DestIterator::column_iterator cd = dupperleft.columnIterator();
1194 ik, ka, kleft, kright, border);
1198 template <
class SrcIterator,
class SrcAccessor,
1199 class DestIterator,
class DestAccessor,
1200 class KernelIterator,
class KernelAccessor>
1203 pair<DestIterator, DestAccessor> dest,
1204 tuple5<KernelIterator, KernelAccessor,
1205 int,
int, BorderTreatmentMode> kernel)
1208 dest.first, dest.second,
1209 kernel.first, kernel.second,
1210 kernel.third, kernel.fourth, kernel.fifth);
1266 template <
class ARITHTYPE>
1307 : iter_(i), base_(i),
1308 count_(count), sum_(count),
1314 throw(PreconditionViolation)
1317 vigra_precondition(count_ == 1 || count_ == sum_,
1318 "Kernel1D::initExplicitly(): "
1319 "Wrong number of init values.");
1344 static value_type one() {
return NumericTraits<value_type>::one(); }
1354 border_treatment_(BORDER_TREATMENT_REFLECT),
1357 kernel_.push_back(norm_);
1363 : kernel_(k.kernel_),
1366 border_treatment_(k.border_treatment_),
1377 border_treatment_(k.borderTreatment()),
1389 border_treatment_ = k.border_treatment_;
1391 kernel_ = k.kernel_;
1414 int size = right_ - left_ + 1;
1415 for(
unsigned int i=0; i<kernel_.
size(); ++i) kernel_[i] = v;
1416 norm_ = (double)size*v;
1418 return InitProxy(kernel_.
begin(),
size, norm_);
1634 this->
initExplicitly(-2, 2) = 0.03134, 0.24, 0.45732, 0.24, 0.03134;
1662 this->
initExplicitly(-2, 2) = 0.04255, 0.241, 0.4329, 0.241, 0.04255;
1690 this->
initExplicitly(-2, 2) = 0.0243, 0.23556, 0.48028, 0.23556, 0.0243;
1725 vigra_precondition(a >= 0.0 && a <= 0.125,
1726 "Kernel1D::initBurtFilter(): 0 <= a <= 0.125 required.");
1964 this->
initExplicitly(-2, 2) = 0.22075, 0.117, -0.6755, 0.117, 0.22075;
2004 vigra_precondition(left <= 0,
2005 "Kernel1D::initExplicitly(): left border must be <= 0.");
2006 vigra_precondition(right >= 0,
2007 "Kernel1D::initExplicitly(): right border must be >= 0.");
2012 kernel_.resize(right - left + 1);
2045 return kernel_[location -
left()];
2050 return kernel_[location -
left()];
2063 int size()
const {
return right_ - left_ + 1; }
2068 {
return border_treatment_; }
2073 { border_treatment_ = new_mode; }
2102 InternalVector kernel_;
2104 BorderTreatmentMode border_treatment_;
2108 template <
class ARITHTYPE>
2110 unsigned int derivativeOrder,
2113 typedef typename NumericTraits<value_type>::RealPromote TmpType;
2117 TmpType sum = NumericTraits<TmpType>::zero();
2119 if(derivativeOrder == 0)
2121 for(; k < kernel_.end(); ++k)
2128 unsigned int faculty = 1;
2129 for(
unsigned int i = 2; i <= derivativeOrder; ++i)
2131 for(
double x = left() + offset; k < kernel_.end(); ++x, ++k)
2133 sum = TmpType(sum + *k * VIGRA_CSTD::pow(-x,
int(derivativeOrder)) / faculty);
2137 vigra_precondition(sum != NumericTraits<value_type>::zero(),
2138 "Kernel1D<ARITHTYPE>::normalize(): "
2139 "Cannot normalize a kernel with sum = 0");
2142 k = kernel_.begin();
2143 for(; k != kernel_.end(); ++k)
2153 template <
class ARITHTYPE>
2158 vigra_precondition(std_dev >= 0.0,
2159 "Kernel1D::initGaussian(): Standard deviation must be >= 0.");
2160 vigra_precondition(windowRatio >= 0.0,
2161 "Kernel1D::initGaussian(): windowRatio must be >= 0.");
2169 if (windowRatio == 0.0)
2170 radius = (int)(3.0 * std_dev + 0.5);
2172 radius = (int)(windowRatio * std_dev + 0.5);
2177 kernel_.erase(kernel_.begin(), kernel_.end());
2178 kernel_.reserve(radius*2+1);
2180 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2182 kernel_.push_back(gauss(x));
2189 kernel_.erase(kernel_.begin(), kernel_.end());
2190 kernel_.push_back(1.0);
2201 border_treatment_ = BORDER_TREATMENT_REFLECT;
2206 template <
class ARITHTYPE>
2210 vigra_precondition(std_dev >= 0.0,
2211 "Kernel1D::initDiscreteGaussian(): Standard deviation must be >= 0.");
2216 int radius = (int)(3.0*std_dev + 0.5);
2220 double f = 2.0 / std_dev / std_dev;
2223 int maxIndex = (int)(2.0 * (radius + 5.0 *
VIGRA_CSTD::sqrt((
double)radius)) + 0.5);
2225 warray[maxIndex] = 0.0;
2226 warray[maxIndex-1] = 1.0;
2228 for(
int i = maxIndex-2; i >= radius; --i)
2230 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2231 if(warray[i] > 1.0e40)
2233 warray[i+1] /= warray[i];
2241 warray[radius+1] = er * warray[radius+1] / warray[radius];
2242 warray[radius] = er;
2244 for(
int i = radius-1; i >= 0; --i)
2246 warray[i] = warray[i+2] + f * (i+1) * warray[i+1];
2250 double scale = norm / (2*er - warray[0]);
2252 initExplicitly(-radius, radius);
2255 for(
int i=0; i<=radius; ++i)
2257 c[i] = c[-i] = warray[i] * scale;
2262 kernel_.erase(kernel_.begin(), kernel_.end());
2263 kernel_.push_back(norm);
2271 border_treatment_ = BORDER_TREATMENT_REFLECT;
2276 template <
class ARITHTYPE>
2283 vigra_precondition(order >= 0,
2284 "Kernel1D::initGaussianDerivative(): Order must be >= 0.");
2288 initGaussian(std_dev, norm, windowRatio);
2292 vigra_precondition(std_dev > 0.0,
2293 "Kernel1D::initGaussianDerivative(): "
2294 "Standard deviation must be > 0.");
2295 vigra_precondition(windowRatio >= 0.0,
2296 "Kernel1D::initGaussianDerivative(): windowRatio must be >= 0.");
2302 if(windowRatio == 0.0)
2303 radius = (int)(3.0 * std_dev + 0.5 * order + 0.5);
2305 radius = (int)(windowRatio * std_dev + 0.5);
2311 kernel_.reserve(radius*2+1);
2316 for(ARITHTYPE x = -(ARITHTYPE)radius; x <= (ARITHTYPE)radius; ++x)
2318 kernel_.push_back(gauss(x));
2319 dc += kernel_[kernel_.size()-1];
2321 dc = ARITHTYPE(dc / (2.0*radius + 1.0));
2327 for(
unsigned int i=0; i < kernel_.size(); ++i)
2337 normalize(norm, order);
2343 border_treatment_ = BORDER_TREATMENT_REFLECT;
2348 template <
class ARITHTYPE>
2353 vigra_precondition(radius > 0,
2354 "Kernel1D::initBinomial(): Radius must be > 0.");
2358 typename InternalVector::iterator x = kernel_.begin() + radius;
2362 for(
int j=radius-1; j>=-radius; --j)
2364 x[j] = 0.5 * x[j+1];
2365 for(
int i=j+1; i<radius; ++i)
2367 x[i] = 0.5 * (x[i] + x[i+1]);
2377 border_treatment_ = BORDER_TREATMENT_REFLECT;
2382 template <
class ARITHTYPE>
2386 vigra_precondition(radius > 0,
2387 "Kernel1D::initAveraging(): Radius must be > 0.");
2390 double scale = 1.0 / (radius * 2 + 1);
2393 kernel_.erase(kernel_.begin(), kernel_.end());
2394 kernel_.reserve(radius*2+1);
2396 for(
int i=0; i<=radius*2+1; ++i)
2398 kernel_.push_back(scale * norm);
2406 border_treatment_ = BORDER_TREATMENT_CLIP;
2411 template <
class ARITHTYPE>
2415 kernel_.erase(kernel_.begin(), kernel_.end());
2418 kernel_.push_back(ARITHTYPE(0.5 * norm));
2419 kernel_.push_back(ARITHTYPE(0.0 * norm));
2420 kernel_.push_back(ARITHTYPE(-0.5 * norm));
2428 border_treatment_ = BORDER_TREATMENT_REFLECT;
2439 template <
class KernelIterator,
class KernelAccessor>
2441 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>
2442 kernel1d(KernelIterator ik, KernelAccessor ka,
2443 int kleft,
int kright, BorderTreatmentMode border)
2446 tuple5<KernelIterator, KernelAccessor, int, int, BorderTreatmentMode>(
2447 ik, ka, kleft, kright, border);
2452 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2453 int, int, BorderTreatmentMode>
2454 kernel1d(Kernel1D<T>
const & k)
2458 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2459 int, int, BorderTreatmentMode>(
2462 k.left(), k.right(),
2463 k.borderTreatment());
2468 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2469 int, int, BorderTreatmentMode>
2470 kernel1d(Kernel1D<T>
const & k, BorderTreatmentMode border)
2474 tuple5<typename Kernel1D<T>::const_iterator,
typename Kernel1D<T>::ConstAccessor,
2475 int, int, BorderTreatmentMode>(
2478 k.left(), k.right(),
2485 #endif // VIGRA_SEPARABLECONVOLUTION_HXX
StandardAccessor< ARITHTYPE > Accessor
Definition: separableconvolution.hxx:1298
InternalVector::const_reference const_reference
Definition: separableconvolution.hxx:1282
void normalize()
Definition: separableconvolution.hxx:2088
StandardConstAccessor< ARITHTYPE > ConstAccessor
Definition: separableconvolution.hxx:1302
value_type norm() const
Definition: separableconvolution.hxx:2077
void initAveraging(int radius, value_type norm)
Definition: separableconvolution.hxx:2383
void initOptimalFirstDerivative5()
Definition: separableconvolution.hxx:1933
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:1267
void initBinomial(int radius, value_type norm)
Definition: separableconvolution.hxx:2350
Kernel1D(Kernel1D const &k)
Definition: separableconvolution.hxx:1362
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2154
InternalVector::value_type value_type
Definition: separableconvolution.hxx:1274
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void initGaussianDerivative(double std_dev, int order)
Definition: separableconvolution.hxx:1525
const_iterator begin() const
Definition: array_vector.hxx:223
InternalVector::iterator iterator
Definition: separableconvolution.hxx:1290
void initGaussian(double std_dev)
Definition: separableconvolution.hxx:1455
void initOptimalSecondDerivative5()
Definition: separableconvolution.hxx:1962
void initDiscreteGaussian(double std_dev)
Definition: separableconvolution.hxx:1483
void initSymmetricDifference()
Definition: separableconvolution.hxx:1879
Kernel1D()
Definition: separableconvolution.hxx:1350
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:2072
void initForwardDifference()
Definition: separableconvolution.hxx:1831
int left() const
Definition: separableconvolution.hxx:2055
Kernel1D & operator=(Kernel1D const &k)
Definition: separableconvolution.hxx:1383
int size() const
Definition: separableconvolution.hxx:2063
void initSecondDifference3()
Definition: separableconvolution.hxx:1901
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
Kernel1D & initExplicitly(int left, int right)
Definition: separableconvolution.hxx:2002
Definition: gaussians.hxx:63
void initDiscreteGaussian(double std_dev, value_type norm)
Definition: separableconvolution.hxx:2207
void initOptimalFirstDerivativeSmoothing3()
Definition: separableconvolution.hxx:1578
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:1683
void initBurtFilter(double a=0.04785)
Definition: separableconvolution.hxx:1723
void initSymmetricGradient()
Definition: separableconvolution.hxx:1808
Accessor accessor()
Definition: separableconvolution.hxx:2099
void initSymmetricGradient(value_type norm)
Definition: separableconvolution.hxx:1798
Kernel1D(Kernel1D< U > const &k)
Definition: separableconvolution.hxx:1373
void initOptimalSecondDerivativeSmoothing5()
Definition: separableconvolution.hxx:1688
InternalVector::reference reference
Definition: separableconvolution.hxx:1278
void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2278
void initOptimalSmoothing5()
Definition: separableconvolution.hxx:1632
Encapsulate read access to the values an iterator points to.
Definition: accessor.hxx:268
~Kernel1D()
Definition: separableconvolution.hxx:1423
BorderTreatmentMode borderTreatment() const
Definition: separableconvolution.hxx:2067
InternalVector::iterator Iterator
Definition: separableconvolution.hxx:1286
int right() const
Definition: separableconvolution.hxx:2059
InternalVector::const_iterator const_iterator
Definition: separableconvolution.hxx:1294
iterator center()
Definition: separableconvolution.hxx:2025
void initOptimalFirstDerivativeSmoothing5()
Definition: separableconvolution.hxx:1660
size_type size() const
Definition: array_vector.hxx:330
void initOptimalSecondDerivativeSmoothing3()
Definition: separableconvolution.hxx:1606
void initAveraging(int radius)
Definition: separableconvolution.hxx:1778
void initBinomial(int radius)
Definition: separableconvolution.hxx:1752
Encapsulate access to the values an iterator points to.
Definition: accessor.hxx:132
InitProxy operator=(value_type const &v)
Definition: separableconvolution.hxx:1412
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:2043
ConstAccessor accessor() const
Definition: separableconvolution.hxx:2095
void initBackwardDifference()
Definition: separableconvolution.hxx:1855
void initOptimalSmoothing3()
Definition: separableconvolution.hxx:1550