[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

nonlineardiffusion.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2002 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_NONLINEARDIFFUSION_HXX
37 #define VIGRA_NONLINEARDIFFUSION_HXX
38 
39 #include <vector>
40 #include "stdimage.hxx"
41 #include "stdimagefunctions.hxx"
42 #include "imageiteratoradapter.hxx"
43 #include "functortraits.hxx"
44 
45 namespace vigra {
46 
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,
52  DestIterator dbegin)
53 {
54  int w = send - sbegin - 1;
55 
56  int i;
57 
58  for(i=0; i<w; ++i)
59  {
60  lower[i] = lower[i] / diag[i];
61 
62  diag[i+1] = diag[i+1] - lower[i] * upper[i];
63  }
64 
65  dbegin[0] = sa(sbegin);
66 
67  for(i=1; i<=w; ++i)
68  {
69  dbegin[i] = sa(sbegin, i) - lower[i-1] * dbegin[i-1];
70  }
71 
72  dbegin[w] = dbegin[w] / diag[w];
73 
74  for(i=w-1; i>=0; --i)
75  {
76  dbegin[i] = (dbegin[i] - upper[i] * dbegin[i+1]) / diag[i];
77  }
78 }
79 
80 
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)
88 {
89  // use traits to determine SumType as to prevent possible overflow
90  typedef typename
91  NumericTraits<typename DestAccessor::value_type>::RealPromote
92  DestType;
93 
94  typedef typename
95  NumericTraits<typename WeightAccessor::value_type>::RealPromote
96  WeightType;
97 
98  // calculate width and height of the image
99  int w = slr.x - sul.x;
100  int h = slr.y - sul.y;
101  int d = (w < h) ? h : w;
102 
103  std::vector<WeightType> lower(d),
104  diag(d),
105  upper(d),
106  res(d);
107 
108  int x,y;
109 
110  WeightType one = NumericTraits<WeightType>::one();
111 
112  // create y iterators
113  SrcIterator ys = sul;
114  WeightIterator yw = wul;
115  DestIterator yd = dul;
116 
117  // x-direction
118  for(y=0; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
119  {
120  typename SrcIterator::row_iterator xs = ys.rowIterator();
121  typename WeightIterator::row_iterator xw = yw.rowIterator();
122  typename DestIterator::row_iterator xd = yd.rowIterator();
123 
124  // fill 3-diag matrix
125 
126  diag[0] = one + timestep * (aw(xw) + aw(xw, 1));
127  for(x=1; x<w-1; ++x)
128  {
129  diag[x] = one + timestep * (2.0 * aw(xw, x) + aw(xw, x+1) + aw(xw, x-1));
130  }
131  diag[w-1] = one + timestep * (aw(xw, w-1) + aw(xw, w-2));
132 
133  for(x=0; x<w-1; ++x)
134  {
135  lower[x] = -timestep * (aw(xw, x) + aw(xw, x+1));
136  upper[x] = lower[x];
137  }
138 
139  internalNonlinearDiffusionDiagonalSolver(xs, xs+w, as,
140  diag.begin(), upper.begin(), lower.begin(), res.begin());
141 
142  for(x=0; x<w; ++x, ++xd)
143  {
144  ad.set(res[x], xd);
145  }
146  }
147 
148  // y-direction
149  ys = sul;
150  yw = wul;
151  yd = dul;
152 
153  for(x=0; x<w; ++x, ++ys.x, ++yd.x, ++yw.x)
154  {
155  typename SrcIterator::column_iterator xs = ys.columnIterator();
156  typename WeightIterator::column_iterator xw = yw.columnIterator();
157  typename DestIterator::column_iterator xd = yd.columnIterator();
158 
159  // fill 3-diag matrix
160 
161  diag[0] = one + timestep * (aw(xw) + aw(xw, 1));
162  for(y=1; y<h-1; ++y)
163  {
164  diag[y] = one + timestep * (2.0 * aw(xw, y) + aw(xw, y+1) + aw(xw, y-1));
165  }
166  diag[h-1] = one + timestep * (aw(xw, h-1) + aw(xw, h-2));
167 
168  for(y=0; y<h-1; ++y)
169  {
170  lower[y] = -timestep * (aw(xw, y) + aw(xw, y+1));
171  upper[y] = lower[y];
172  }
173 
174  internalNonlinearDiffusionDiagonalSolver(xs, xs+h, as,
175  diag.begin(), upper.begin(), lower.begin(), res.begin());
176 
177  for(y=0; y<h; ++y, ++xd)
178  {
179  ad.set(0.5 * (ad(xd) + res[y]), xd);
180  }
181  }
182 }
183 
184 /** \addtogroup NonLinearDiffusion Non-linear Diffusion
185 
186  Perform edge-preserving smoothing.
187 */
188 //@{
189 
190 /********************************************************/
191 /* */
192 /* nonlinearDiffusion */
193 /* */
194 /********************************************************/
195 
196 /** \brief Perform edge-preserving smoothing at the given scale.
197 
198  The algorithm solves the non-linear diffusion equation
199 
200  \f[
201  \frac{\partial}{\partial t} u =
202  \frac{\partial}{\partial x}
203  \left( g(|\nabla u|)
204  \frac{\partial}{\partial x} u
205  \right)
206  \f]
207 
208  where <em> t</em> is the time, <b> x</b> is the location vector,
209  <em> u(</em><b> x</b><em> , t)</em> is the smoothed image at time <em> t</em>, and
210  <em> g(.)</em> is the location dependent diffusivity. At time zero, the image
211  <em> u(</em><b> x</b><em> , 0)</em> is simply the original image. The time is
212  proportional to the square of the scale parameter: \f$t = s^2\f$.
213  The diffusion
214  equation is solved iteratively according
215  to the Additive Operator Splitting Scheme (AOS) from
216 
217 
218  J. Weickert: <em>"Recursive Separable Schemes for Nonlinear Diffusion
219  Filters"</em>,
220  in: B. ter Haar Romeny, L. Florack, J. Koenderingk, M. Viergever (eds.):
221  1st Intl. Conf. on Scale-Space Theory in Computer Vision 1997,
222  Springer LNCS 1252
223 
224  <TT>DiffusivityFunctor</TT> implements the gradient dependent local diffusivity.
225  It is passed
226  as an argument to \ref gradientBasedTransform(). The return value must be
227  between 0 and 1 and determines the weight a pixel gets when
228  its neighbors are smoothed. Weickert recommends the use of the diffusivity
229  implemented by class \ref DiffusivityFunctor. It's also possible to use
230  other functors, for example one that always returns 1, in which case
231  we obtain the solution to the linear diffusion equation, i.e.
232  Gaussian convolution.
233 
234  The source value type must be a
235  linear space with internal addition, scalar multiplication, and
236  NumericTraits defined. The value_type of the DiffusivityFunctor must be the
237  scalar field over wich the source value type's linear space is defined.
238 
239  In addition to <TT>nonlinearDiffusion()</TT>, there is an algorithm
240  <TT>nonlinearDiffusionExplicit()</TT> which implements the Explicit Scheme
241  described in the above article. Both algorithms have the same interface,
242  but the explicit scheme gives slightly more accurate approximations of
243  the diffusion process at the cost of much slower processing.
244 
245  <b> Declarations:</b>
246 
247  pass arguments explicitly:
248  \code
249  namespace vigra {
250  template <class SrcIterator, class SrcAccessor,
251  class DestIterator, class DestAccessor,
252  class DiffusivityFunctor>
253  void nonlinearDiffusion(SrcIterator sul, SrcIterator slr, SrcAccessor as,
254  DestIterator dul, DestAccessor ad,
255  DiffusivityFunctor const & weight, double scale);
256  }
257  \endcode
258 
259 
260  use argument objects in conjunction with \ref ArgumentObjectFactories :
261  \code
262  namespace vigra {
263  template <class SrcIterator, class SrcAccessor,
264  class DestIterator, class DestAccessor,
265  class DiffusivityFunctor>
266  void nonlinearDiffusion(
267  triple<SrcIterator, SrcIterator, SrcAccessor> src,
268  pair<DestIterator, DestAccessor> dest,
269  DiffusivityFunctor const & weight, double scale);
270  }
271  \endcode
272 
273  <b> Usage:</b>
274 
275  <b>\#include</b> <vigra/nonlineardiffusion.hxx>
276 
277 
278  \code
279  FImage src(w,h), dest(w,h);
280  float edge_threshold, scale;
281  ...
282 
283  nonlinearDiffusion(srcImageRange(src), destImage(dest),
284  DiffusivityFunctor<float>(edge_threshold), scale);
285  \endcode
286 
287  <b> Required Interface:</b>
288 
289  <ul>
290 
291  <li> <TT>SrcIterator</TT> and <TT>DestIterator</TT> are models of ImageIterator
292  <li> <TT>SrcAccessor</TT> and <TT>DestAccessor</TT> are models of StandardAccessor
293  <li> <TT>SrcAccessor::value_type</TT> is a linear space
294  <li> <TT>DiffusivityFunctor</TT> conforms to the requirements of
295  \ref gradientBasedTransform(). Its range is between 0 and 1.
296  <li> <TT>DiffusivityFunctor::value_type</TT> is an algebraic field
297 
298  </ul>
299 
300  <b> Precondition:</b>
301 
302  <TT>scale > 0</TT>
303 */
304 doxygen_overloaded_function(template <...> void nonlinearDiffusion)
305 
306 template <class SrcIterator, class SrcAccessor,
307  class DestIterator, class DestAccessor,
308  class DiffusivityFunc>
309 void nonlinearDiffusion(SrcIterator sul, SrcIterator slr, SrcAccessor as,
310  DestIterator dul, DestAccessor ad,
311  DiffusivityFunc const & weight, double scale)
312 {
313  vigra_precondition(scale > 0.0, "nonlinearDiffusion(): scale must be > 0");
314 
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;
319 
320  Size2D size(slr.x - sul.x, slr.y - sul.y);
321 
322  typedef typename
323  NumericTraits<typename SrcAccessor::value_type>::RealPromote
324  TmpType;
325  typedef typename DiffusivityFunc::value_type WeightType;
326 
327  BasicImage<TmpType> smooth1(size);
328  BasicImage<TmpType> smooth2(size);
329 
330  BasicImage<WeightType> weights(size);
331 
332  typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
333  s2 = smooth2.upperLeft();
334  typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
335 
336  typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
337  typename BasicImage<WeightType>::Accessor wa = weights.accessor();
338 
339  gradientBasedTransform(sul, slr, as, wi, wa, weight);
340 
341  internalNonlinearDiffusionAOSStep(sul, slr, as, wi, wa, s1, a, rest_time);
342 
343  for(int i = 0; i < number_of_steps; ++i)
344  {
345  gradientBasedTransform(s1, s1+size, a, wi, wa, weight);
346 
347  internalNonlinearDiffusionAOSStep(s1, s1+size, a, wi, wa, s2, a, time_step);
348 
349  std::swap(s1, s2);
350  }
351 
352  copyImage(s1, s1+size, a, dul, ad);
353 }
354 
355 template <class SrcIterator, class SrcAccessor,
356  class DestIterator, class DestAccessor,
357  class DiffusivityFunc>
358 inline
359 void nonlinearDiffusion(
360  triple<SrcIterator, SrcIterator, SrcAccessor> src,
361  pair<DestIterator, DestAccessor> dest,
362  DiffusivityFunc const & weight, double scale)
363 {
364  nonlinearDiffusion(src.first, src.second, src.third,
365  dest.first, dest.second,
366  weight, scale);
367 }
368 
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,
376  double time_step)
377 {
378  // use traits to determine SumType as to prevent possible overflow
379  typedef typename
380  NumericTraits<typename SrcAccessor::value_type>::RealPromote
381  SumType;
382 
383  typedef typename
384  NumericTraits<typename WeightAccessor::value_type>::RealPromote
385  WeightType;
386 
387  // calculate width and height of the image
388  int w = slr.x - sul.x;
389  int h = slr.y - sul.y;
390 
391  int x,y;
392 
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);
397 
398  WeightType gt, gb, gl, gr, g0;
399  WeightType one = NumericTraits<WeightType>::one();
400  SumType sum;
401 
402  time_step /= 2.0;
403 
404  // create y iterators
405  SrcIterator ys = sul;
406  WeightIterator yw = wul;
407  DestIterator yd = dul;
408 
409  SrcIterator xs = ys;
410  WeightIterator xw = yw;
411  DestIterator xd = yd;
412 
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;
418 
419  sum = g0 * as(xs);
420  sum += gt * as(xs, bottom);
421  sum += gb * as(xs, bottom);
422  sum += gl * as(xs, right);
423  sum += gr * as(xs, right);
424 
425  ad.set(sum, xd);
426 
427  for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
428  {
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;
434 
435  sum = g0 * as(xs);
436  sum += gt * as(xs, bottom);
437  sum += gb * as(xs, bottom);
438  sum += gl * as(xs, left);
439  sum += gr * as(xs, right);
440 
441  ad.set(sum, xd);
442  }
443 
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;
449 
450  sum = g0 * as(xs);
451  sum += gt * as(xs, bottom);
452  sum += gb * as(xs, bottom);
453  sum += gl * as(xs, left);
454  sum += gr * as(xs, left);
455 
456  ad.set(sum, xd);
457 
458  for(y=2, ++ys.y, ++yd.y, ++yw.y; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
459  {
460  xs = ys;
461  xd = yd;
462  xw = yw;
463 
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;
469 
470  sum = g0 * as(xs);
471  sum += gt * as(xs, top);
472  sum += gb * as(xs, bottom);
473  sum += gl * as(xs, right);
474  sum += gr * as(xs, right);
475 
476  ad.set(sum, xd);
477 
478  for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
479  {
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;
485 
486  sum = g0 * as(xs);
487  sum += gt * as(xs, top);
488  sum += gb * as(xs, bottom);
489  sum += gl * as(xs, left);
490  sum += gr * as(xs, right);
491 
492  ad.set(sum, xd);
493  }
494 
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;
500 
501  sum = g0 * as(xs);
502  sum += gt * as(xs, top);
503  sum += gb * as(xs, bottom);
504  sum += gl * as(xs, left);
505  sum += gr * as(xs, left);
506 
507  ad.set(sum, xd);
508  }
509 
510  xs = ys;
511  xd = yd;
512  xw = yw;
513 
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;
519 
520  sum = g0 * as(xs);
521  sum += gt * as(xs, top);
522  sum += gb * as(xs, top);
523  sum += gl * as(xs, right);
524  sum += gr * as(xs, right);
525 
526  ad.set(sum, xd);
527 
528  for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
529  {
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;
535 
536  sum = g0 * as(xs);
537  sum += gt * as(xs, top);
538  sum += gb * as(xs, top);
539  sum += gl * as(xs, left);
540  sum += gr * as(xs, right);
541 
542  ad.set(sum, xd);
543  }
544 
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;
550 
551  sum = g0 * as(xs);
552  sum += gt * as(xs, top);
553  sum += gb * as(xs, top);
554  sum += gl * as(xs, left);
555  sum += gr * as(xs, left);
556 
557  ad.set(sum, xd);
558 }
559 
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)
566 {
567  vigra_precondition(scale > 0.0, "nonlinearDiffusionExplicit(): scale must be > 0");
568 
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;
573 
574  Size2D size(slr.x - sul.x, slr.y - sul.y);
575 
576  typedef typename
577  NumericTraits<typename SrcAccessor::value_type>::RealPromote
578  TmpType;
579  typedef typename DiffusivityFunc::value_type WeightType;
580 
581  BasicImage<TmpType> smooth1(size);
582  BasicImage<TmpType> smooth2(size);
583 
584  BasicImage<WeightType> weights(size);
585 
586  typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
587  s2 = smooth2.upperLeft();
588  typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
589 
590  typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
591  typename BasicImage<WeightType>::Accessor wa = weights.accessor();
592 
593  gradientBasedTransform(sul, slr, as, wi, wa, weight);
594 
595  internalNonlinearDiffusionExplicitStep(sul, slr, as, wi, wa, s1, a, rest_time);
596 
597  for(int i = 0; i < number_of_steps; ++i)
598  {
599  gradientBasedTransform(s1, s1+size, a, wi, wa, weight);
600 
601  internalNonlinearDiffusionExplicitStep(s1, s1+size, a, wi, wa, s2, a, time_step);
602 
603  swap(s1, s2);
604  }
605 
606  copyImage(s1, s1+size, a, dul, ad);
607 }
608 
609 template <class SrcIterator, class SrcAccessor,
610  class DestIterator, class DestAccessor,
611  class DiffusivityFunc>
612 inline
613 void nonlinearDiffusionExplicit(
614  triple<SrcIterator, SrcIterator, SrcAccessor> src,
615  pair<DestIterator, DestAccessor> dest,
616  DiffusivityFunc const & weight, double scale)
617 {
618  nonlinearDiffusionExplicit(src.first, src.second, src.third,
619  dest.first, dest.second,
620  weight, scale);
621 }
622 
623 /********************************************************/
624 /* */
625 /* DiffusivityFunctor */
626 /* */
627 /********************************************************/
628 
629 /** \brief Diffusivity functor for non-linear diffusion.
630 
631  This functor implements the diffusivity recommended by Weickert:
632 
633  \f[
634  g(|\nabla u|) = 1 -
635  \exp{\left(\frac{-3.315}{(|\nabla u| / thresh)^4}\right)}
636  \f]
637 
638 
639  where <TT>thresh</TT> is a threshold that determines whether a specific gradient
640  magnitude is interpreted as a significant edge (above the threshold)
641  or as noise. It is meant to be used with \ref nonlinearDiffusion().
642 */
643 template <class Value>
645 {
646  public:
647  /** the functors first argument type (must be an algebraic field with <TT>exp()</TT> defined).
648  However, the functor also works with RGBValue<first_argument_type> (this hack is
649  necessary since Microsoft C++ doesn't support partial specialization).
650  */
651  typedef Value first_argument_type;
652 
653  /** the functors second argument type (same as first).
654  However, the functor also works with RGBValue<second_argument_type> (this hack is
655  necessary since Microsoft C++ doesn't support partial specialization).
656  */
657  typedef Value second_argument_type;
658 
659  /** the functors result type
660  */
661  typedef typename NumericTraits<Value>::RealPromote result_type;
662 
663  /** \deprecated use first_argument_type, second_argument_type, result_type
664  */
665  typedef Value value_type;
666 
667  /** init functor with given edge threshold
668  */
669  DiffusivityFunctor(Value const & thresh)
670  : weight_(thresh*thresh),
671  one_(NumericTraits<result_type>::one()),
672  zero_(NumericTraits<result_type>::zero())
673  {}
674 
675  /** calculate diffusivity from scalar arguments
676  */
678  {
679  Value mag = (gx*gx + gy*gy) / weight_;
680 
681  return (mag == zero_) ? one_ : one_ - VIGRA_CSTD::exp(-3.315 / mag / mag);
682  }
683 
684  /** calculate diffusivity from RGB arguments
685  */
687  {
688  result_type mag = (gx.red()*gx.red() +
689  gx.green()*gx.green() +
690  gx.blue()*gx.blue() +
691  gy.red()*gy.red() +
692  gy.green()*gy.green() +
693  gy.blue()*gy.blue()) / weight_;
694 
695  return (mag == zero_) ? one_ : one_ - VIGRA_CSTD::exp(-3.315 / mag / mag);
696  }
697 
698  result_type weight_;
699  result_type one_;
700  result_type zero_;
701 };
702 
703 template <class ValueType>
704 class FunctorTraits<DiffusivityFunctor<ValueType> >
705 : public FunctorTraitsBase<DiffusivityFunctor<ValueType> >
706 {
707  public:
708  typedef VigraTrueType isBinaryFunctor;
709 };
710 
711 
712 //@}
713 
714 } // namespace vigra
715 
716 #endif /* VIGRA_NONLINEARDIFFUSION_HXX */

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (Wed Sep 26 2012)