36 #ifndef VIGRA_NONLINEARDIFFUSION_HXX
37 #define VIGRA_NONLINEARDIFFUSION_HXX
40 #include "stdimage.hxx"
41 #include "stdimagefunctions.hxx"
42 #include "imageiteratoradapter.hxx"
43 #include "functortraits.hxx"
44 #include "multi_shape.hxx"
48 template <
class SrcIterator,
class SrcAccessor,
49 class CoeffIterator,
class DestIterator>
50 void internalNonlinearDiffusionDiagonalSolver(
51 SrcIterator sbegin, SrcIterator send, SrcAccessor sa,
52 CoeffIterator diag, CoeffIterator upper, CoeffIterator lower,
55 int w = send - sbegin - 1;
61 lower[i] = lower[i] / diag[i];
63 diag[i+1] = diag[i+1] - lower[i] * upper[i];
66 dbegin[0] = sa(sbegin);
70 dbegin[i] = sa(sbegin, i) - lower[i-1] * dbegin[i-1];
73 dbegin[w] = dbegin[w] / diag[w];
77 dbegin[i] = (dbegin[i] - upper[i] * dbegin[i+1]) / diag[i];
82 template <
class SrcIterator,
class SrcAccessor,
83 class WeightIterator,
class WeightAccessor,
84 class DestIterator,
class DestAccessor>
85 void internalNonlinearDiffusionAOSStep(
86 SrcIterator sul, SrcIterator slr, SrcAccessor as,
87 WeightIterator wul, WeightAccessor aw,
88 DestIterator dul, DestAccessor ad,
double timestep)
92 NumericTraits<typename WeightAccessor::value_type>::RealPromote
96 int w = slr.x - sul.x;
97 int h = slr.y - sul.y;
98 int d = (w < h) ? h : w;
100 std::vector<WeightType> lower(d),
107 WeightType one = NumericTraits<WeightType>::one();
110 SrcIterator ys = sul;
111 WeightIterator yw = wul;
112 DestIterator yd = dul;
115 for(y=0; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
117 typename SrcIterator::row_iterator xs = ys.rowIterator();
118 typename WeightIterator::row_iterator xw = yw.rowIterator();
119 typename DestIterator::row_iterator xd = yd.rowIterator();
123 diag[0] = one + timestep * (aw(xw) + aw(xw, 1));
126 diag[x] = one + timestep * (2.0 * aw(xw, x) + aw(xw, x+1) + aw(xw, x-1));
128 diag[w-1] = one + timestep * (aw(xw, w-1) + aw(xw, w-2));
132 lower[x] = -timestep * (aw(xw, x) + aw(xw, x+1));
136 internalNonlinearDiffusionDiagonalSolver(xs, xs+w, as,
137 diag.begin(), upper.begin(), lower.begin(), res.begin());
139 for(x=0; x<w; ++x, ++xd)
150 for(x=0; x<w; ++x, ++ys.x, ++yd.x, ++yw.x)
152 typename SrcIterator::column_iterator xs = ys.columnIterator();
153 typename WeightIterator::column_iterator xw = yw.columnIterator();
154 typename DestIterator::column_iterator xd = yd.columnIterator();
158 diag[0] = one + timestep * (aw(xw) + aw(xw, 1));
161 diag[y] = one + timestep * (2.0 * aw(xw, y) + aw(xw, y+1) + aw(xw, y-1));
163 diag[h-1] = one + timestep * (aw(xw, h-1) + aw(xw, h-2));
167 lower[y] = -timestep * (aw(xw, y) + aw(xw, y+1));
171 internalNonlinearDiffusionDiagonalSolver(xs, xs+h, as,
172 diag.begin(), upper.begin(), lower.begin(), res.begin());
174 for(y=0; y<h; ++y, ++xd)
176 ad.set(0.5 * (ad(xd) + res[y]), xd);
331 template <
class SrcIterator,
class SrcAccessor,
332 class DestIterator,
class DestAccessor,
333 class DiffusivityFunc>
335 DestIterator dul, DestAccessor ad,
336 DiffusivityFunc
const & weight,
double scale)
338 vigra_precondition(scale > 0.0,
"nonlinearDiffusion(): scale must be > 0");
340 double total_time = scale*scale/2.0;
341 const double time_step = 5.0;
342 int number_of_steps = (int)(total_time / time_step);
343 double rest_time = total_time - time_step * number_of_steps;
345 Size2D size(slr.x - sul.x, slr.y - sul.y);
348 NumericTraits<typename SrcAccessor::value_type>::RealPromote
350 typedef typename DiffusivityFunc::value_type WeightType;
352 BasicImage<TmpType> smooth1(size);
353 BasicImage<TmpType> smooth2(size);
355 BasicImage<WeightType> weights(size);
357 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
358 s2 = smooth2.upperLeft();
359 typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
361 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
362 typename BasicImage<WeightType>::Accessor wa = weights.accessor();
366 internalNonlinearDiffusionAOSStep(sul, slr, as, wi, wa, s1, a, rest_time);
368 for(
int i = 0; i < number_of_steps; ++i)
372 internalNonlinearDiffusionAOSStep(s1, s1+size, a, wi, wa, s2, a, time_step);
380 template <
class SrcIterator,
class SrcAccessor,
381 class DestIterator,
class DestAccessor,
382 class DiffusivityFunc>
385 pair<DestIterator, DestAccessor> dest,
386 DiffusivityFunc
const & weight,
double scale)
389 dest.first, dest.second,
393 template <
class T1,
class S1,
395 class DiffusivityFunc>
398 MultiArrayView<2, T2, S2> dest,
399 DiffusivityFunc
const & weight,
double scale)
401 vigra_precondition(src.shape() == dest.shape(),
402 "nonlinearDiffusion(): shape mismatch between input and output.");
410 template <
class SrcIterator,
class SrcAccessor,
411 class WeightIterator,
class WeightAccessor,
412 class DestIterator,
class DestAccessor>
413 void internalNonlinearDiffusionExplicitStep(
414 SrcIterator sul, SrcIterator slr, SrcAccessor as,
415 WeightIterator wul, WeightAccessor aw,
416 DestIterator dul, DestAccessor ad,
421 NumericTraits<typename SrcAccessor::value_type>::RealPromote
425 NumericTraits<typename WeightAccessor::value_type>::RealPromote
429 int w = slr.x - sul.x;
430 int h = slr.y - sul.y;
434 const Diff2D left(-1, 0);
435 const Diff2D right(1, 0);
436 const Diff2D top(0, -1);
437 const Diff2D bottom(0, 1);
439 WeightType gt, gb, gl, gr, g0;
440 WeightType one = NumericTraits<WeightType>::one();
446 SrcIterator ys = sul;
447 WeightIterator yw = wul;
448 DestIterator yd = dul;
451 WeightIterator xw = yw;
452 DestIterator xd = yd;
454 gt = (aw(xw) + aw(xw, bottom)) * time_step;
455 gb = (aw(xw) + aw(xw, bottom)) * time_step;
456 gl = (aw(xw) + aw(xw, right)) * time_step;
457 gr = (aw(xw) + aw(xw, right)) * time_step;
458 g0 = one - gt - gb - gr - gl;
461 sum += gt * as(xs, bottom);
462 sum += gb * as(xs, bottom);
463 sum += gl * as(xs, right);
464 sum += gr * as(xs, right);
468 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
470 gt = (aw(xw) + aw(xw, bottom)) * time_step;
471 gb = (aw(xw) + aw(xw, bottom)) * time_step;
472 gl = (aw(xw) + aw(xw, left)) * time_step;
473 gr = (aw(xw) + aw(xw, right)) * time_step;
474 g0 = one - gt - gb - gr - gl;
477 sum += gt * as(xs, bottom);
478 sum += gb * as(xs, bottom);
479 sum += gl * as(xs, left);
480 sum += gr * as(xs, right);
485 gt = (aw(xw) + aw(xw, bottom)) * time_step;
486 gb = (aw(xw) + aw(xw, bottom)) * time_step;
487 gl = (aw(xw) + aw(xw, left)) * time_step;
488 gr = (aw(xw) + aw(xw, left)) * time_step;
489 g0 = one - gt - gb - gr - gl;
492 sum += gt * as(xs, bottom);
493 sum += gb * as(xs, bottom);
494 sum += gl * as(xs, left);
495 sum += gr * as(xs, left);
499 for(y=2, ++ys.y, ++yd.y, ++yw.y; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
505 gt = (aw(xw) + aw(xw, top)) * time_step;
506 gb = (aw(xw) + aw(xw, bottom)) * time_step;
507 gl = (aw(xw) + aw(xw, right)) * time_step;
508 gr = (aw(xw) + aw(xw, right)) * time_step;
509 g0 = one - gt - gb - gr - gl;
512 sum += gt * as(xs, top);
513 sum += gb * as(xs, bottom);
514 sum += gl * as(xs, right);
515 sum += gr * as(xs, right);
519 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
521 gt = (aw(xw) + aw(xw, top)) * time_step;
522 gb = (aw(xw) + aw(xw, bottom)) * time_step;
523 gl = (aw(xw) + aw(xw, left)) * time_step;
524 gr = (aw(xw) + aw(xw, right)) * time_step;
525 g0 = one - gt - gb - gr - gl;
528 sum += gt * as(xs, top);
529 sum += gb * as(xs, bottom);
530 sum += gl * as(xs, left);
531 sum += gr * as(xs, right);
536 gt = (aw(xw) + aw(xw, top)) * time_step;
537 gb = (aw(xw) + aw(xw, bottom)) * time_step;
538 gl = (aw(xw) + aw(xw, left)) * time_step;
539 gr = (aw(xw) + aw(xw, left)) * time_step;
540 g0 = one - gt - gb - gr - gl;
543 sum += gt * as(xs, top);
544 sum += gb * as(xs, bottom);
545 sum += gl * as(xs, left);
546 sum += gr * as(xs, left);
555 gt = (aw(xw) + aw(xw, top)) * time_step;
556 gb = (aw(xw) + aw(xw, top)) * time_step;
557 gl = (aw(xw) + aw(xw, right)) * time_step;
558 gr = (aw(xw) + aw(xw, right)) * time_step;
559 g0 = one - gt - gb - gr - gl;
562 sum += gt * as(xs, top);
563 sum += gb * as(xs, top);
564 sum += gl * as(xs, right);
565 sum += gr * as(xs, right);
569 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
571 gt = (aw(xw) + aw(xw, top)) * time_step;
572 gb = (aw(xw) + aw(xw, top)) * time_step;
573 gl = (aw(xw) + aw(xw, left)) * time_step;
574 gr = (aw(xw) + aw(xw, right)) * time_step;
575 g0 = one - gt - gb - gr - gl;
578 sum += gt * as(xs, top);
579 sum += gb * as(xs, top);
580 sum += gl * as(xs, left);
581 sum += gr * as(xs, right);
586 gt = (aw(xw) + aw(xw, top)) * time_step;
587 gb = (aw(xw) + aw(xw, top)) * time_step;
588 gl = (aw(xw) + aw(xw, left)) * time_step;
589 gr = (aw(xw) + aw(xw, left)) * time_step;
590 g0 = one - gt - gb - gr - gl;
593 sum += gt * as(xs, top);
594 sum += gb * as(xs, top);
595 sum += gl * as(xs, left);
596 sum += gr * as(xs, left);
607 template <
class SrcIterator,
class SrcAccessor,
608 class DestIterator,
class DestAccessor,
609 class DiffusivityFunc>
611 DestIterator dul, DestAccessor ad,
612 DiffusivityFunc
const & weight,
double scale)
614 vigra_precondition(scale > 0.0,
"nonlinearDiffusionExplicit(): scale must be > 0");
616 double total_time = scale*scale/2.0;
617 const double time_step = 0.25;
618 int number_of_steps = total_time / time_step;
619 double rest_time = total_time - time_step * number_of_steps;
621 Size2D size(slr.x - sul.x, slr.y - sul.y);
624 NumericTraits<typename SrcAccessor::value_type>::RealPromote
626 typedef typename DiffusivityFunc::value_type WeightType;
628 BasicImage<TmpType> smooth1(size);
629 BasicImage<TmpType> smooth2(size);
631 BasicImage<WeightType> weights(size);
633 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
634 s2 = smooth2.upperLeft();
635 typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
637 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
638 typename BasicImage<WeightType>::Accessor wa = weights.accessor();
642 internalNonlinearDiffusionExplicitStep(sul, slr, as, wi, wa, s1, a, rest_time);
644 for(
int i = 0; i < number_of_steps; ++i)
648 internalNonlinearDiffusionExplicitStep(s1, s1+size, a, wi, wa, s2, a, time_step);
656 template <
class SrcIterator,
class SrcAccessor,
657 class DestIterator,
class DestAccessor,
658 class DiffusivityFunc>
661 pair<DestIterator, DestAccessor> dest,
662 DiffusivityFunc
const & weight,
double scale)
665 dest.first, dest.second,
669 template <
class T1,
class S1,
671 class DiffusivityFunc>
674 MultiArrayView<2, T2, S2> dest,
675 DiffusivityFunc
const & weight,
double scale)
677 vigra_precondition(src.shape() == dest.shape(),
678 "nonlinearDiffusionExplicit(): shape mismatch between input and output.");
704 template <
class Value>
731 : weight_(thresh*thresh),
740 Value mag = (gx*gx + gy*gy) / weight_;
742 return (mag == zero_) ? one_ : one_ -
VIGRA_CSTD::exp(-3.315 / mag / mag);
756 return (mag == zero_) ? one_ : one_ -
VIGRA_CSTD::exp(-3.315 / mag / mag);
764 template <
class ValueType>
765 class FunctorTraits<DiffusivityFunctor<ValueType> >
766 :
public FunctorTraitsBase<DiffusivityFunctor<ValueType> >
769 typedef VigraTrueType isBinaryFunctor;
Value first_argument_type
Definition: nonlineardiffusion.hxx:712
Diffusivity functor for non-linear diffusion.
Definition: nonlineardiffusion.hxx:705
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
value_type & red()
Definition: rgbvalue.hxx:278
Value value_type
Definition: nonlineardiffusion.hxx:726
Value second_argument_type
Definition: nonlineardiffusion.hxx:718
value_type & blue()
Definition: rgbvalue.hxx:286
DiffusivityFunctor(Value const &thresh)
Definition: nonlineardiffusion.hxx:730
value_type & green()
Definition: rgbvalue.hxx:282
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void copyImage(...)
Copy source image into destination image.
result_type operator()(RGBValue< Value > const &gx, RGBValue< Value > const &gy) const
Definition: nonlineardiffusion.hxx:747
void nonlinearDiffusionExplicit(...)
Perform edge-preserving smoothing at the given scale using an explicit scheme.
void nonlinearDiffusion(...)
Perform edge-preserving smoothing at the given scale.
result_type operator()(first_argument_type const &gx, second_argument_type const &gy) const
Definition: nonlineardiffusion.hxx:738
Class for a single RGB value.
Definition: accessor.hxx:938
NumericTraits< Value >::RealPromote result_type
Definition: nonlineardiffusion.hxx:722