37 #ifndef VIGRA_STDCONVOLUTION_HXX
38 #define VIGRA_STDCONVOLUTION_HXX
41 #include "stdimage.hxx"
42 #include "bordertreatment.hxx"
43 #include "separableconvolution.hxx"
44 #include "utilities.hxx"
45 #include "sized_int.hxx"
179 template <
class SrcIterator,
class SrcAccessor,
180 class DestIterator,
class DestAccessor,
181 class KernelIterator,
class KernelAccessor>
182 void convolveImage(SrcIterator src_ul, SrcIterator src_lr, SrcAccessor src_acc,
183 DestIterator dest_ul, DestAccessor dest_acc,
184 KernelIterator ki, KernelAccessor ak,
187 vigra_precondition((border == BORDER_TREATMENT_CLIP ||
188 border == BORDER_TREATMENT_AVOID ||
189 border == BORDER_TREATMENT_REFLECT ||
190 border == BORDER_TREATMENT_REPEAT ||
191 border == BORDER_TREATMENT_WRAP),
193 " Border treatment must be one of follow treatments:\n"
194 " - BORDER_TREATMENT_CLIP\n"
195 " - BORDER_TREATMENT_AVOID\n"
196 " - BORDER_TREATMENT_REFLECT\n"
197 " - BORDER_TREATMENT_REPEAT\n"
198 " - BORDER_TREATMENT_WRAP\n");
200 vigra_precondition(kul.
x <= 0 && kul.
y <= 0,
201 "convolveImage(): coordinates of "
202 "kernel's upper left must be <= 0.");
203 vigra_precondition(klr.
x >= 0 && klr.
y >= 0,
204 "convolveImage(): coordinates of "
205 "kernel's lower right must be >= 0.");
209 PromoteTraits<
typename SrcAccessor::value_type,
210 typename KernelAccessor::value_type>::Promote SumType;
212 NumericTraits<typename KernelAccessor::value_type>::RealPromote KernelSumType;
213 typedef typename DestAccessor::value_type DestType;
216 int w = src_lr.x - src_ul.x;
217 int h = src_lr.y - src_ul.y;
220 int kernel_width = klr.
x - kul.
x + 1;
221 int kernel_height = klr.
y - kul.
y + 1;
224 "convolveImage(): kernel larger than image.");
226 KernelSumType
norm = NumericTraits<KernelSumType>::zero();
227 if(border == BORDER_TREATMENT_CLIP)
230 KernelIterator yk = ki + klr;
233 for(
int y = 0; y < kernel_height; ++y, --yk.y)
235 KernelIterator xk = yk;
236 for(
int x = 0; x < kernel_width; ++x, --xk.x)
241 vigra_precondition(norm != NumericTraits<KernelSumType>::zero(),
242 "convolveImage(): Cannot use BORDER_TREATMENT_CLIP with a DC-free kernel");
245 DestIterator yd = dest_ul;
246 SrcIterator ys = src_ul;
247 SrcIterator send = src_lr;
250 for(
int y=0; y<h; ++y, ++ys.y, ++yd.y)
256 for(
int x=0; x < w; ++x, ++xs.x, ++xd.x)
259 SumType
sum = NumericTraits<SumType>::zero();
260 KernelIterator ykernel = ki + klr;
262 if(x >= klr.x && y >= klr.y && x < w + kul.
x && y < h + kul.
y)
265 SrcIterator yys = xs - klr;
266 SrcIterator yyend = xs - kul;
268 for(; yys.y <= yyend.y; ++yys.y, --ykernel.y)
270 typename SrcIterator::row_iterator xxs = yys.
rowIterator();
271 typename SrcIterator::row_iterator xxe = xxs + kernel_width;
272 typename KernelIterator::row_iterator xkernel= ykernel.rowIterator();
274 for(; xxs < xxe; ++xxs, --xkernel)
276 sum += ak(xkernel) * src_acc(xxs);
280 else if(border == BORDER_TREATMENT_REPEAT)
283 for(
int yk = klr.y; yk >= kul.
y; --yk, --ykernel.y)
286 typename KernelIterator::row_iterator xkernel = ykernel.rowIterator();
288 for(
int xk = klr.x; xk >= kul.
x; --xk, --xkernel)
291 sum += ak(xkernel) * src_acc(src_ul, diff);
295 else if(border == BORDER_TREATMENT_REFLECT)
298 for(
int yk = klr.y; yk >= kul.
y; --yk , --ykernel.y)
300 diff.
y =
abs(y - yk);
302 diff.
y = 2*h - 2 - diff.
y;
303 typename KernelIterator::row_iterator xkernel = ykernel.rowIterator();
305 for(
int xk = klr.x; xk >= kul.
x; --xk, --xkernel)
307 diff.
x =
abs(x - xk);
309 diff.
x = 2*w - 2 - diff.
x;
310 sum += ak(xkernel) * src_acc(src_ul, diff);
314 else if(border == BORDER_TREATMENT_WRAP)
317 for(
int yk = klr.y; yk >= kul.
y; --yk, --ykernel.y)
319 diff.
y = (y - yk + h) % h;
320 typename KernelIterator::row_iterator xkernel = ykernel.rowIterator();
322 for(
int xk = klr.x; xk >= kul.
x; --xk, --xkernel)
324 diff.
x = (x - xk + w) % w;
325 sum += ak(xkernel) * src_acc(src_ul, diff);
329 else if(border == BORDER_TREATMENT_CLIP)
331 KernelSumType ksum = NumericTraits<KernelSumType>::zero();
333 for(
int yk = klr.y; yk >= kul.
y; --yk, --ykernel.y)
336 if(diff.
y < 0 || diff.
y >= h)
338 typename KernelIterator::row_iterator xkernel = ykernel.rowIterator();
340 for(
int xk = klr.x; xk >= kul.
x; --xk, --xkernel)
343 if(diff.
x < 0 || diff.
x >= w)
346 sum += ak(xkernel) * src_acc(src_ul, diff);
352 else if(border == BORDER_TREATMENT_AVOID)
358 dest_acc.set(detail::RequiresExplicitCast<DestType>::cast(sum), xd);
363 template <
class SrcIterator,
class SrcAccessor,
364 class DestIterator,
class DestAccessor,
365 class KernelIterator,
class KernelAccessor>
368 triple<SrcIterator, SrcIterator, SrcAccessor> src,
369 pair<DestIterator, DestAccessor> dest,
370 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
371 BorderTreatmentMode> kernel)
374 dest.first, dest.second,
375 kernel.first, kernel.second, kernel.third,
376 kernel.fourth, kernel.fifth);
521 template <
class SrcIterator,
class SrcAccessor,
522 class DestIterator,
class DestAccessor,
523 class MaskIterator,
class MaskAccessor,
524 class KernelIterator,
class KernelAccessor>
527 MaskIterator
mul, MaskAccessor am,
528 DestIterator dest_ul, DestAccessor dest_acc,
529 KernelIterator ki, KernelAccessor ak,
530 Diff2D kul, Diff2D klr, BorderTreatmentMode border)
532 vigra_precondition((border == BORDER_TREATMENT_CLIP ||
533 border == BORDER_TREATMENT_AVOID),
534 "normalizedConvolveImage(): "
535 "Border treatment must be BORDER_TREATMENT_CLIP or BORDER_TREATMENT_AVOID.");
537 vigra_precondition(kul.x <= 0 && kul.y <= 0,
538 "normalizedConvolveImage(): left borders must be <= 0.");
539 vigra_precondition(klr.x >= 0 && klr.y >= 0,
540 "normalizedConvolveImage(): right borders must be >= 0.");
544 NumericTraits<typename SrcAccessor::value_type>::RealPromote SumType;
546 NumericTraits<typename KernelAccessor::value_type>::RealPromote KSumType;
548 NumericTraits<typename DestAccessor::value_type> DestTraits;
551 int w = src_lr.x - src_ul.x;
552 int h = src_lr.y - src_ul.y;
553 int kernel_width = klr.x - kul.x + 1;
554 int kernel_height = klr.y - kul.y + 1;
557 int ystart = (border == BORDER_TREATMENT_AVOID) ? klr.y : 0;
558 int yend = (border == BORDER_TREATMENT_AVOID) ? h+kul.y : h;
559 int xstart = (border == BORDER_TREATMENT_AVOID) ? klr.x : 0;
560 int xend = (border == BORDER_TREATMENT_AVOID) ? w+kul.x : w;
563 DestIterator yd = dest_ul + Diff2D(xstart, ystart);
564 SrcIterator ys = src_ul + Diff2D(xstart, ystart);
565 MaskIterator ym = mul + Diff2D(xstart, ystart);
567 KSumType
norm = ak(ki);
569 KernelIterator yk = ki + klr;
570 for(yy=0; yy<kernel_height; ++yy, --yk.y)
572 KernelIterator xk = yk;
574 for(xx=0; xx<kernel_width; ++xx, --xk.x)
582 for(y=ystart; y < yend; ++y, ++ys.y, ++yd.y, ++ym.y)
589 for(x=xstart; x < xend; ++x, ++xs.x, ++xd.x, ++xm.x)
594 y0 = (y<klr.y) ? -y : -klr.y;
595 y1 = (h-y-1<-kul.y) ? h-y-1 : -kul.y;
596 x0 = (x<klr.x) ? -x : -klr.x;
597 x1 = (w-x-1<-kul.x) ? w-x-1 : -kul.x;
601 SumType
sum = NumericTraits<SumType>::zero();
602 KSumType ksum = NumericTraits<KSumType>::zero();
604 SrcIterator yys = xs + Diff2D(x0, y0);
605 MaskIterator yym = xm + Diff2D(x0, y0);
606 KernelIterator yk = ki - Diff2D(x0, y0);
608 int xx, kernel_width, kernel_height;
609 kernel_width = x1 - x0 + 1;
610 kernel_height = y1 - y0 + 1;
611 for(yy=0; yy<kernel_height; ++yy, ++yys.y, --yk.y, ++yym.y)
613 typename SrcIterator::row_iterator xxs = yys.rowIterator();
614 typename SrcIterator::row_iterator xxend = xxs + kernel_width;
615 typename MaskIterator::row_iterator xxm = yym.rowIterator();
616 typename KernelIterator::row_iterator xk = yk.rowIterator();
618 for(xx=0; xxs < xxend; ++xxs, --xk, ++xxm)
620 if(!am(xxm))
continue;
624 sum = detail::RequiresExplicitCast<SumType>::cast(ak(xk) * src_acc(xxs));
630 sum = detail::RequiresExplicitCast<SumType>::cast(sum + ak(xk) * src_acc(xxs));
636 if(ksum != NumericTraits<KSumType>::zero())
638 dest_acc.set(DestTraits::fromRealPromote(
639 detail::RequiresExplicitCast<SumType>::cast((norm / ksum) * sum)), xd);
646 template <
class SrcIterator,
class SrcAccessor,
647 class DestIterator,
class DestAccessor,
648 class MaskIterator,
class MaskAccessor,
649 class KernelIterator,
class KernelAccessor>
652 triple<SrcIterator, SrcIterator, SrcAccessor> src,
653 pair<MaskIterator, MaskAccessor> mask,
654 pair<DestIterator, DestAccessor> dest,
655 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
656 BorderTreatmentMode> kernel)
659 mask.first, mask.second,
660 dest.first, dest.second,
661 kernel.first, kernel.second, kernel.third,
662 kernel.fourth, kernel.fifth);
705 template <
class SrcIterator,
class SrcAccessor,
706 class DestIterator,
class DestAccessor,
707 class MaskIterator,
class MaskAccessor,
708 class KernelIterator,
class KernelAccessor>
711 MaskIterator mul, MaskAccessor am,
712 DestIterator dest_ul, DestAccessor dest_acc,
713 KernelIterator ki, KernelAccessor ak,
714 Diff2D kul, Diff2D klr, BorderTreatmentMode border)
719 ki, ak, kul, klr, border);
722 template <
class SrcIterator,
class SrcAccessor,
723 class DestIterator,
class DestAccessor,
724 class MaskIterator,
class MaskAccessor,
725 class KernelIterator,
class KernelAccessor>
728 triple<SrcIterator, SrcIterator, SrcAccessor> src,
729 pair<MaskIterator, MaskAccessor> mask,
730 pair<DestIterator, DestAccessor> dest,
731 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D,
732 BorderTreatmentMode> kernel)
735 mask.first, mask.second,
736 dest.first, dest.second,
737 kernel.first, kernel.second, kernel.third,
738 kernel.fourth, kernel.fifth);
797 template <
class ARITHTYPE>
827 : iter_(i), base_(i),
828 count_(count), sum_(count),
834 vigra_precondition(count_ == 1 || count_ == sum_,
835 "Kernel2D::initExplicitly(): "
836 "Too few init values.");
841 if(count_ == sum_) norm_ = *iter_;
844 vigra_precondition(count_ > 0,
845 "Kernel2D::initExplicitly(): "
846 "Too many init values.");
861 static value_type one() {
return NumericTraits<value_type>::one(); }
868 : kernel_(1, 1, one()),
872 border_treatment_(BORDER_TREATMENT_REFLECT)
878 : kernel_(k.kernel_),
882 border_treatment_(k.border_treatment_)
895 border_treatment_ = k.border_treatment_;
921 int size = (right_.
x - left_.
x + 1) *
922 (right_.
y - left_.
y + 1);
924 norm_ = (double)size*v;
926 return InitProxy(kernel_.
begin(), size, norm_);
952 int w = right_.
x - left_.
x + 1;
953 int h = right_.
y - left_.
y + 1;
961 KIter kiy = ky.
center() + left_.
y;
964 for(
int y=left_.
y; y<=right_.
y; ++y, ++kiy, ++iy.y)
966 KIter kix = kx.
center() + left_.
x;
968 for(
int x=left_.
x; x<=right_.
x; ++x, ++kix, ++ix.x)
970 *ix = ka(kix) * ka(kiy);
998 template <
class KernelIterator>
1000 KernelIterator kycenter,
int yleft,
int yright)
1002 vigra_precondition(xleft <= 0 && yleft <= 0,
1003 "Kernel2D::initSeparable(): left borders must be <= 0.");
1004 vigra_precondition(xright >= 0 && yright >= 0,
1005 "Kernel2D::initSeparable(): right borders must be >= 0.");
1007 left_ =
Point2D(xleft, yleft);
1008 right_ =
Point2D(xright, yright);
1010 int w = right_.
x - left_.
x + 1;
1011 int h = right_.
y - left_.
y + 1;
1014 KernelIterator kiy = kycenter + left_.
y;
1017 for(
int y=left_.
y; y<=right_.
y; ++y, ++kiy, ++iy.y)
1019 KernelIterator kix = kxcenter + left_.
x;
1021 for(
int x=left_.
x; x<=right_.
x; ++x, ++kix, ++ix.x)
1032 for(; i!= iend; ++i)
1051 initGaussian(std_dev, NumericTraits<value_type>::one());
1076 vigra_precondition(radius > 0,
1077 "Kernel2D::initDisk(): radius must be > 0.");
1079 left_ =
Point2D(-radius, -radius);
1080 right_ =
Point2D(radius, radius);
1081 int w = right_.
x - left_.
x + 1;
1082 int h = right_.
y - left_.
y + 1;
1084 norm_ = NumericTraits<value_type>::one();
1086 kernel_ = NumericTraits<value_type>::zero();
1090 double r2 = (double)radius*radius;
1093 for(i=0; i<= radius; ++i)
1095 double r = (double) i - 0.5;
1097 for(
int j=-w; j<=w; ++j)
1099 k(j, i) = NumericTraits<value_type>::one();
1100 k(j, -i) = NumericTraits<value_type>::one();
1101 count += (i != 0) ? 2.0 : 1.0;
1105 count = 1.0 / count;
1107 for(
int y=-radius; y<=radius; ++y)
1109 for(
int x=-radius; x<=radius; ++x)
1111 k(x,y) = count * k(x,y);
1156 vigra_precondition(upperleft.
x <= 0 && upperleft.
y <= 0,
1157 "Kernel2D::initExplicitly(): left borders must be <= 0.");
1158 vigra_precondition(lowerright.
x >= 0 && lowerright.
y >= 0,
1159 "Kernel2D::initExplicitly(): right borders must be >= 0.");
1164 int w = right_.
x - left_.
x + 1;
1165 int h = right_.
y - left_.
y + 1;
1181 int width()
const {
return right_.
x - left_.
x + 1; }
1198 {
return kernel_[
Diff2D(x,y) - left_]; }
1203 {
return kernel_[
Diff2D(x,y) - left_]; }
1208 {
return kernel_[d - left_]; }
1213 {
return kernel_[d - left_]; }
1246 typename NumericTraits<value_type>::RealPromote sum = *i;
1249 for(; i!= iend; ++i)
1255 i = kernel_.
begin();
1256 for(; i != iend; ++i)
1274 {
return border_treatment_; }
1282 vigra_precondition((new_mode == BORDER_TREATMENT_CLIP ||
1283 new_mode == BORDER_TREATMENT_AVOID ||
1284 new_mode == BORDER_TREATMENT_REFLECT ||
1285 new_mode == BORDER_TREATMENT_REPEAT ||
1286 new_mode == BORDER_TREATMENT_WRAP),
1287 "convolveImage():\n"
1288 " Border treatment must be one of follow treatments:\n"
1289 " - BORDER_TREATMENT_CLIP\n"
1290 " - BORDER_TREATMENT_AVOID\n"
1291 " - BORDER_TREATMENT_REFLECT\n"
1292 " - BORDER_TREATMENT_REPEAT\n"
1293 " - BORDER_TREATMENT_WRAP\n");
1295 border_treatment_ = new_mode;
1303 BorderTreatmentMode border_treatment_;
1314 template <
class KernelIterator,
class KernelAccessor>
1316 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode>
1317 kernel2d(KernelIterator ik, KernelAccessor ak, Diff2D kul, Diff2D klr,
1318 BorderTreatmentMode border)
1322 tuple5<KernelIterator, KernelAccessor, Diff2D, Diff2D, BorderTreatmentMode> (
1323 ik, ak, kul, klr, border);
1328 tuple5<typename Kernel2D<T>::ConstIterator,
1329 typename Kernel2D<T>::ConstAccessor,
1330 Diff2D, Diff2D, BorderTreatmentMode>
1331 kernel2d(Kernel2D<T>
const & k)
1335 tuple5<typename Kernel2D<T>::ConstIterator,
1336 typename Kernel2D<T>::ConstAccessor,
1337 Diff2D, Diff2D, BorderTreatmentMode>(
1340 k.upperLeft(), k.lowerRight(),
1341 k.borderTreatment());
1346 tuple5<typename Kernel2D<T>::ConstIterator,
1347 typename Kernel2D<T>::ConstAccessor,
1348 Diff2D, Diff2D, BorderTreatmentMode>
1349 kernel2d(Kernel2D<T>
const & k, BorderTreatmentMode border)
1353 tuple5<typename Kernel2D<T>::ConstIterator,
1354 typename Kernel2D<T>::ConstAccessor,
1355 Diff2D, Diff2D, BorderTreatmentMode>(
1358 k.upperLeft(), k.lowerRight(),
1365 #endif // VIGRA_STDCONVOLUTION_HXX