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