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"
47 template <
class SrcIterator,
class SrcAccessor,
48 class CoeffIterator,
class DestIterator>
49 void internalNonlinearDiffusionDiagonalSolver(
50 SrcIterator sbegin, SrcIterator send, SrcAccessor sa,
51 CoeffIterator diag, CoeffIterator upper, CoeffIterator lower,
54 int w = send - sbegin - 1;
60 lower[i] = lower[i] / diag[i];
62 diag[i+1] = diag[i+1] - lower[i] * upper[i];
65 dbegin[0] = sa(sbegin);
69 dbegin[i] = sa(sbegin, i) - lower[i-1] * dbegin[i-1];
72 dbegin[w] = dbegin[w] / diag[w];
76 dbegin[i] = (dbegin[i] - upper[i] * dbegin[i+1]) / diag[i];
81 template <
class SrcIterator,
class SrcAccessor,
82 class WeightIterator,
class WeightAccessor,
83 class DestIterator,
class DestAccessor>
84 void internalNonlinearDiffusionAOSStep(
85 SrcIterator sul, SrcIterator slr, SrcAccessor as,
86 WeightIterator wul, WeightAccessor aw,
87 DestIterator dul, DestAccessor ad,
double timestep)
91 NumericTraits<typename DestAccessor::value_type>::RealPromote
95 NumericTraits<typename WeightAccessor::value_type>::RealPromote
99 int w = slr.x - sul.x;
100 int h = slr.y - sul.y;
101 int d = (w < h) ? h : w;
103 std::vector<WeightType> lower(d),
110 WeightType one = NumericTraits<WeightType>::one();
113 SrcIterator ys = sul;
114 WeightIterator yw = wul;
115 DestIterator yd = dul;
118 for(y=0; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
120 typename SrcIterator::row_iterator xs = ys.rowIterator();
121 typename WeightIterator::row_iterator xw = yw.rowIterator();
122 typename DestIterator::row_iterator xd = yd.rowIterator();
126 diag[0] = one + timestep * (aw(xw) + aw(xw, 1));
129 diag[x] = one + timestep * (2.0 * aw(xw, x) + aw(xw, x+1) + aw(xw, x-1));
131 diag[w-1] = one + timestep * (aw(xw, w-1) + aw(xw, w-2));
135 lower[x] = -timestep * (aw(xw, x) + aw(xw, x+1));
139 internalNonlinearDiffusionDiagonalSolver(xs, xs+w, as,
140 diag.begin(), upper.begin(), lower.begin(), res.begin());
142 for(x=0; x<w; ++x, ++xd)
153 for(x=0; x<w; ++x, ++ys.x, ++yd.x, ++yw.x)
155 typename SrcIterator::column_iterator xs = ys.columnIterator();
156 typename WeightIterator::column_iterator xw = yw.columnIterator();
157 typename DestIterator::column_iterator xd = yd.columnIterator();
161 diag[0] = one + timestep * (aw(xw) + aw(xw, 1));
164 diag[y] = one + timestep * (2.0 * aw(xw, y) + aw(xw, y+1) + aw(xw, y-1));
166 diag[h-1] = one + timestep * (aw(xw, h-1) + aw(xw, h-2));
170 lower[y] = -timestep * (aw(xw, y) + aw(xw, y+1));
174 internalNonlinearDiffusionDiagonalSolver(xs, xs+h, as,
175 diag.begin(), upper.begin(), lower.begin(), res.begin());
177 for(y=0; y<h; ++y, ++xd)
179 ad.set(0.5 * (ad(xd) + res[y]), xd);
306 template <
class SrcIterator,
class SrcAccessor,
307 class DestIterator,
class DestAccessor,
308 class DiffusivityFunc>
310 DestIterator dul, DestAccessor ad,
311 DiffusivityFunc
const & weight,
double scale)
313 vigra_precondition(scale > 0.0,
"nonlinearDiffusion(): scale must be > 0");
315 double total_time = scale*scale/2.0;
316 static const double time_step = 5.0;
317 int number_of_steps = (int)(total_time / time_step);
318 double rest_time = total_time - time_step * number_of_steps;
320 Size2D size(slr.x - sul.x, slr.y - sul.y);
323 NumericTraits<typename SrcAccessor::value_type>::RealPromote
325 typedef typename DiffusivityFunc::value_type WeightType;
327 BasicImage<TmpType> smooth1(size);
328 BasicImage<TmpType> smooth2(size);
330 BasicImage<WeightType> weights(size);
332 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
333 s2 = smooth2.upperLeft();
334 typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
336 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
337 typename BasicImage<WeightType>::Accessor wa = weights.accessor();
341 internalNonlinearDiffusionAOSStep(sul, slr, as, wi, wa, s1, a, rest_time);
343 for(
int i = 0; i < number_of_steps; ++i)
347 internalNonlinearDiffusionAOSStep(s1, s1+size, a, wi, wa, s2, a, time_step);
355 template <
class SrcIterator,
class SrcAccessor,
356 class DestIterator,
class DestAccessor,
357 class DiffusivityFunc>
360 triple<SrcIterator, SrcIterator, SrcAccessor> src,
361 pair<DestIterator, DestAccessor> dest,
362 DiffusivityFunc
const & weight,
double scale)
365 dest.first, dest.second,
369 template <
class SrcIterator,
class SrcAccessor,
370 class WeightIterator,
class WeightAccessor,
371 class DestIterator,
class DestAccessor>
372 void internalNonlinearDiffusionExplicitStep(
373 SrcIterator sul, SrcIterator slr, SrcAccessor as,
374 WeightIterator wul, WeightAccessor aw,
375 DestIterator dul, DestAccessor ad,
380 NumericTraits<typename SrcAccessor::value_type>::RealPromote
384 NumericTraits<typename WeightAccessor::value_type>::RealPromote
388 int w = slr.x - sul.x;
389 int h = slr.y - sul.y;
393 static const Diff2D left(-1, 0);
394 static const Diff2D right(1, 0);
395 static const Diff2D top(0, -1);
396 static const Diff2D bottom(0, 1);
398 WeightType gt, gb, gl, gr, g0;
399 WeightType one = NumericTraits<WeightType>::one();
405 SrcIterator ys = sul;
406 WeightIterator yw = wul;
407 DestIterator yd = dul;
410 WeightIterator xw = yw;
411 DestIterator xd = yd;
413 gt = (aw(xw) + aw(xw, bottom)) * time_step;
414 gb = (aw(xw) + aw(xw, bottom)) * time_step;
415 gl = (aw(xw) + aw(xw, right)) * time_step;
416 gr = (aw(xw) + aw(xw, right)) * time_step;
417 g0 = one - gt - gb - gr - gl;
420 sum += gt * as(xs, bottom);
421 sum += gb * as(xs, bottom);
422 sum += gl * as(xs, right);
423 sum += gr * as(xs, right);
427 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
429 gt = (aw(xw) + aw(xw, bottom)) * time_step;
430 gb = (aw(xw) + aw(xw, bottom)) * time_step;
431 gl = (aw(xw) + aw(xw, left)) * time_step;
432 gr = (aw(xw) + aw(xw, right)) * time_step;
433 g0 = one - gt - gb - gr - gl;
436 sum += gt * as(xs, bottom);
437 sum += gb * as(xs, bottom);
438 sum += gl * as(xs, left);
439 sum += gr * as(xs, right);
444 gt = (aw(xw) + aw(xw, bottom)) * time_step;
445 gb = (aw(xw) + aw(xw, bottom)) * time_step;
446 gl = (aw(xw) + aw(xw, left)) * time_step;
447 gr = (aw(xw) + aw(xw, left)) * time_step;
448 g0 = one - gt - gb - gr - gl;
451 sum += gt * as(xs, bottom);
452 sum += gb * as(xs, bottom);
453 sum += gl * as(xs, left);
454 sum += gr * as(xs, left);
458 for(y=2, ++ys.y, ++yd.y, ++yw.y; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
464 gt = (aw(xw) + aw(xw, top)) * time_step;
465 gb = (aw(xw) + aw(xw, bottom)) * time_step;
466 gl = (aw(xw) + aw(xw, right)) * time_step;
467 gr = (aw(xw) + aw(xw, right)) * time_step;
468 g0 = one - gt - gb - gr - gl;
471 sum += gt * as(xs, top);
472 sum += gb * as(xs, bottom);
473 sum += gl * as(xs, right);
474 sum += gr * as(xs, right);
478 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
480 gt = (aw(xw) + aw(xw, top)) * time_step;
481 gb = (aw(xw) + aw(xw, bottom)) * time_step;
482 gl = (aw(xw) + aw(xw, left)) * time_step;
483 gr = (aw(xw) + aw(xw, right)) * time_step;
484 g0 = one - gt - gb - gr - gl;
487 sum += gt * as(xs, top);
488 sum += gb * as(xs, bottom);
489 sum += gl * as(xs, left);
490 sum += gr * as(xs, right);
495 gt = (aw(xw) + aw(xw, top)) * time_step;
496 gb = (aw(xw) + aw(xw, bottom)) * time_step;
497 gl = (aw(xw) + aw(xw, left)) * time_step;
498 gr = (aw(xw) + aw(xw, left)) * time_step;
499 g0 = one - gt - gb - gr - gl;
502 sum += gt * as(xs, top);
503 sum += gb * as(xs, bottom);
504 sum += gl * as(xs, left);
505 sum += gr * as(xs, left);
514 gt = (aw(xw) + aw(xw, top)) * time_step;
515 gb = (aw(xw) + aw(xw, top)) * time_step;
516 gl = (aw(xw) + aw(xw, right)) * time_step;
517 gr = (aw(xw) + aw(xw, right)) * time_step;
518 g0 = one - gt - gb - gr - gl;
521 sum += gt * as(xs, top);
522 sum += gb * as(xs, top);
523 sum += gl * as(xs, right);
524 sum += gr * as(xs, right);
528 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
530 gt = (aw(xw) + aw(xw, top)) * time_step;
531 gb = (aw(xw) + aw(xw, top)) * time_step;
532 gl = (aw(xw) + aw(xw, left)) * time_step;
533 gr = (aw(xw) + aw(xw, right)) * time_step;
534 g0 = one - gt - gb - gr - gl;
537 sum += gt * as(xs, top);
538 sum += gb * as(xs, top);
539 sum += gl * as(xs, left);
540 sum += gr * as(xs, right);
545 gt = (aw(xw) + aw(xw, top)) * time_step;
546 gb = (aw(xw) + aw(xw, top)) * time_step;
547 gl = (aw(xw) + aw(xw, left)) * time_step;
548 gr = (aw(xw) + aw(xw, left)) * time_step;
549 g0 = one - gt - gb - gr - gl;
552 sum += gt * as(xs, top);
553 sum += gb * as(xs, top);
554 sum += gl * as(xs, left);
555 sum += gr * as(xs, left);
560 template <
class SrcIterator,
class SrcAccessor,
561 class DestIterator,
class DestAccessor,
562 class DiffusivityFunc>
563 void nonlinearDiffusionExplicit(SrcIterator sul, SrcIterator slr, SrcAccessor as,
564 DestIterator dul, DestAccessor ad,
565 DiffusivityFunc
const & weight,
double scale)
567 vigra_precondition(scale > 0.0,
"nonlinearDiffusionExplicit(): scale must be > 0");
569 double total_time = scale*scale/2.0;
570 static const double time_step = 0.25;
571 int number_of_steps = total_time / time_step;
572 double rest_time = total_time - time_step * number_of_steps;
574 Size2D size(slr.x - sul.x, slr.y - sul.y);
577 NumericTraits<typename SrcAccessor::value_type>::RealPromote
579 typedef typename DiffusivityFunc::value_type WeightType;
581 BasicImage<TmpType> smooth1(size);
582 BasicImage<TmpType> smooth2(size);
584 BasicImage<WeightType> weights(size);
586 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
587 s2 = smooth2.upperLeft();
588 typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
590 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
591 typename BasicImage<WeightType>::Accessor wa = weights.accessor();
595 internalNonlinearDiffusionExplicitStep(sul, slr, as, wi, wa, s1, a, rest_time);
597 for(
int i = 0; i < number_of_steps; ++i)
601 internalNonlinearDiffusionExplicitStep(s1, s1+size, a, wi, wa, s2, a, time_step);
609 template <
class SrcIterator,
class SrcAccessor,
610 class DestIterator,
class DestAccessor,
611 class DiffusivityFunc>
613 void nonlinearDiffusionExplicit(
614 triple<SrcIterator, SrcIterator, SrcAccessor> src,
615 pair<DestIterator, DestAccessor> dest,
616 DiffusivityFunc
const & weight,
double scale)
618 nonlinearDiffusionExplicit(src.first, src.second, src.third,
619 dest.first, dest.second,
643 template <
class Value>
670 : weight_(thresh*thresh),
679 Value mag = (gx*gx + gy*gy) / weight_;
681 return (mag == zero_) ? one_ : one_ -
VIGRA_CSTD::exp(-3.315 / mag / mag);
695 return (mag == zero_) ? one_ : one_ -
VIGRA_CSTD::exp(-3.315 / mag / mag);
703 template <
class ValueType>
704 class FunctorTraits<DiffusivityFunctor<ValueType> >
705 :
public FunctorTraitsBase<DiffusivityFunctor<ValueType> >
708 typedef VigraTrueType isBinaryFunctor;
Value first_argument_type
Definition: nonlineardiffusion.hxx:651
Diffusivity functor for non-linear diffusion.
Definition: nonlineardiffusion.hxx:644
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
value_type & red()
Definition: rgbvalue.hxx:279
Value value_type
Definition: nonlineardiffusion.hxx:665
Value second_argument_type
Definition: nonlineardiffusion.hxx:657
value_type & blue()
Definition: rgbvalue.hxx:287
DiffusivityFunctor(Value const &thresh)
Definition: nonlineardiffusion.hxx:669
value_type & green()
Definition: rgbvalue.hxx:283
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:1683
void copyImage(...)
Copy source image into destination image.
result_type operator()(RGBValue< Value > const &gx, RGBValue< Value > const &gy) const
Definition: nonlineardiffusion.hxx:686
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:677
Class for a single RGB value.
Definition: accessor.hxx:937
NumericTraits< Value >::RealPromote result_type
Definition: nonlineardiffusion.hxx:661