[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/nonlineardiffusion.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.2.0, Aug 07 2003 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 #ifndef VIGRA_NONLINEARDIFFUSION_HXX 00024 #define VIGRA_NONLINEARDIFFUSION_HXX 00025 00026 #include <vector> 00027 #include "vigra/stdimage.hxx" 00028 #include "vigra/stdimagefunctions.hxx" 00029 #include "vigra/imageiteratoradapter.hxx" 00030 00031 namespace vigra { 00032 00033 template <class SrcIterator, class SrcAccessor, 00034 class CoeffIterator, class DestIterator> 00035 void internalNonlinearDiffusionDiagonalSolver( 00036 SrcIterator sbegin, SrcIterator send, SrcAccessor sa, 00037 CoeffIterator diag, CoeffIterator upper, CoeffIterator lower, 00038 DestIterator dbegin) 00039 { 00040 int w = send - sbegin - 1; 00041 00042 int i; 00043 00044 for(i=0; i<w; ++i) 00045 { 00046 lower[i] = lower[i] / diag[i]; 00047 00048 diag[i+1] = diag[i+1] - lower[i] * upper[i]; 00049 } 00050 00051 dbegin[0] = sa(sbegin); 00052 00053 for(i=1; i<=w; ++i) 00054 { 00055 dbegin[i] = sa(sbegin, i) - lower[i-1] * dbegin[i-1]; 00056 } 00057 00058 dbegin[w] = dbegin[w] / diag[w]; 00059 00060 for(i=w-1; i>=0; --i) 00061 { 00062 dbegin[i] = (dbegin[i] - upper[i] * dbegin[i+1]) / diag[i]; 00063 } 00064 } 00065 00066 00067 template <class SrcIterator, class SrcAccessor, 00068 class WeightIterator, class WeightAccessor, 00069 class DestIterator, class DestAccessor> 00070 void internalNonlinearDiffusionAOSStep( 00071 SrcIterator sul, SrcIterator slr, SrcAccessor as, 00072 WeightIterator wul, WeightAccessor aw, 00073 DestIterator dul, DestAccessor ad, double timestep) 00074 { 00075 // use traits to determine SumType as to prevent possible overflow 00076 typedef typename 00077 NumericTraits<typename DestAccessor::value_type>::RealPromote 00078 DestType; 00079 00080 typedef typename 00081 NumericTraits<typename WeightAccessor::value_type>::RealPromote 00082 WeightType; 00083 00084 // calculate width and height of the image 00085 int w = slr.x - sul.x; 00086 int h = slr.y - sul.y; 00087 int d = (w < h) ? h : w; 00088 00089 std::vector<WeightType> lower(d), 00090 diag(d), 00091 upper(d); 00092 00093 std::vector<DestType> res(d); 00094 00095 int x,y; 00096 00097 WeightType one = NumericTraits<WeightType>::one(); 00098 00099 // create y iterators 00100 SrcIterator ys = sul; 00101 WeightIterator yw = wul; 00102 DestIterator yd = dul; 00103 00104 // x-direction 00105 for(y=0; y<h; ++y, ++ys.y, ++yd.y, ++yw.y) 00106 { 00107 typename SrcIterator::row_iterator xs = ys.rowIterator(); 00108 typename WeightIterator::row_iterator xw = yw.rowIterator(); 00109 typename DestIterator::row_iterator xd = yd.rowIterator(); 00110 00111 // fill 3-diag matrix 00112 00113 diag[0] = one + timestep * (aw(xw) + aw(xw, 1)); 00114 for(x=1; x<w-1; ++x) 00115 { 00116 diag[x] = one + timestep * (2.0 * aw(xw, x) + aw(xw, x+1) + aw(xw, x-1)); 00117 } 00118 diag[w-1] = one + timestep * (aw(xw, w-1) + aw(xw, w-2)); 00119 00120 for(x=0; x<w-1; ++x) 00121 { 00122 lower[x] = -timestep * (aw(xw, x) + aw(xw, x+1)); 00123 upper[x] = lower[x]; 00124 } 00125 00126 internalNonlinearDiffusionDiagonalSolver(xs, xs+w, as, 00127 diag.begin(), upper.begin(), lower.begin(), res.begin()); 00128 00129 for(x=0; x<w; ++x, ++xd) 00130 { 00131 ad.set(res[x], xd); 00132 } 00133 } 00134 00135 // y-direction 00136 ys = sul; 00137 yw = wul; 00138 yd = dul; 00139 00140 for(x=0; x<w; ++x, ++ys.x, ++yd.x, ++yw.x) 00141 { 00142 typename SrcIterator::column_iterator xs = ys.columnIterator(); 00143 typename WeightIterator::column_iterator xw = yw.columnIterator(); 00144 typename DestIterator::column_iterator xd = yd.columnIterator(); 00145 00146 // fill 3-diag matrix 00147 00148 diag[0] = one + timestep * (aw(xw) + aw(xw, 1)); 00149 for(y=1; y<h-1; ++y) 00150 { 00151 diag[y] = one + timestep * (2.0 * aw(xw, y) + aw(xw, y+1) + aw(xw, y-1)); 00152 } 00153 diag[h-1] = one + timestep * (aw(xw, h-1) + aw(xw, h-2)); 00154 00155 for(y=0; y<h-1; ++y) 00156 { 00157 lower[y] = -timestep * (aw(xw, y) + aw(xw, y+1)); 00158 upper[y] = lower[y]; 00159 } 00160 00161 internalNonlinearDiffusionDiagonalSolver(xs, xs+h, as, 00162 diag.begin(), upper.begin(), lower.begin(), res.begin()); 00163 00164 for(y=0; y<h; ++y, ++xd) 00165 { 00166 ad.set(0.5 * (ad(xd) + res[y]), xd); 00167 } 00168 } 00169 } 00170 00171 /** \addtogroup NonLinearDiffusion Non-linear Diffusion 00172 00173 Perform edge-preserving smoothing. 00174 */ 00175 //@{ 00176 00177 /********************************************************/ 00178 /* */ 00179 /* nonlinearDiffusion */ 00180 /* */ 00181 /********************************************************/ 00182 00183 /** \brief Perform edge-preserving smoothing at the given scale. 00184 00185 The algorithm solves the non-linear diffusion equation 00186 00187 \f[ 00188 \frac{\partial}{\partial t} u = 00189 \frac{\partial}{\partial x} 00190 \left( g(|\nabla u|) 00191 \frac{\partial}{\partial x} u 00192 \right) 00193 \f] 00194 00195 where <em> t</em> is the time, <b> x</b> is the location vector, 00196 <em> u(</em><b> x</b><em> , t)</em> is the smoothed image at time <em> t</em>, and 00197 <em> g(.)</em> is the location dependent diffusivity. At time zero, the image 00198 <em> u(</em><b> x</b><em> , 0)</em> is simply the original image. The time is 00199 propotional to the square of the scale parameter: \f$t = s^2\f$. 00200 The diffusion 00201 equation is solved iteratively according 00202 to the Additive Operator Splitting Scheme (AOS) from 00203 00204 00205 J. Weickert: <em>"Recursive Separable Schemes for Nonlinear Diffusion 00206 Filters"</em>, 00207 in: B. ter Haar Romeny, L. Florack, J. Koenderingk, M. Viergever (eds.): 00208 1st Intl. Conf. on Scale-Space Theory in Computer Vision 1997, 00209 Springer LNCS 1252 00210 00211 <TT>DiffusivityFunctor</TT> implements the gradient dependent local diffusivity. 00212 It is passed 00213 as an argument to \ref gradientBasedTransform(). The return value must be 00214 between 0 and 1 and determines the weight a pixel gets when 00215 its neighbors are smoothed. Weickert recommends the use of the diffusivity 00216 implemented by class \ref DiffusivityFunctor. It's also possible to use 00217 other functors, for example one that always returns 1, in which case 00218 we obtain the solution to the linear diffusion equation, i.e. 00219 Gaussian convolution. 00220 00221 The source value type must be a 00222 linear space with internal addition, scalar multiplication, and 00223 NumericTraits defined. The value_type of the DiffusivityFunctor must be the 00224 scalar field over wich the source value type's linear space is defined. 00225 00226 In addition to <TT>nonlinearDiffusion()</TT>, there is an algorithm 00227 <TT>nonlinearDiffusionExplicit()</TT> which implements the Explicit Scheme 00228 described in the above article. Both algorithms have the same interface, 00229 but the explicit scheme gives slightly more accurate approximations of 00230 the diffusion process at the cost of much slower processing. 00231 00232 <b> Declarations:</b> 00233 00234 pass arguments explicitly: 00235 \code 00236 namespace vigra { 00237 template <class SrcIterator, class SrcAccessor, 00238 class DestIterator, class DestAccessor, 00239 class DiffusivityFunctor> 00240 void nonlinearDiffusion(SrcIterator sul, SrcIterator slr, SrcAccessor as, 00241 DestIterator dul, DestAccessor ad, 00242 DiffusivityFunctor const & weight, double scale); 00243 } 00244 \endcode 00245 00246 00247 use argument objects in conjuction with \ref ArgumentObjectFactories: 00248 \code 00249 namespace vigra { 00250 template <class SrcIterator, class SrcAccessor, 00251 class DestIterator, class DestAccessor, 00252 class DiffusivityFunctor> 00253 void nonlinearDiffusion( 00254 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00255 pair<DestIterator, DestAccessor> dest, 00256 DiffusivityFunctor const & weight, double scale); 00257 } 00258 \endcode 00259 00260 <b> Usage:</b> 00261 00262 <b>\#include</b> "<a href="nonlineardiffusion_8hxx-source.html">vigra/nonlineardiffusion.hxx</a>" 00263 00264 00265 \code 00266 FImage src(w,h), dest(w,h); 00267 float edge_threshold, scale; 00268 ... 00269 00270 nonlinearDiffusion(srcImageRange(src), destImage(dest), 00271 DiffusivityFunctor<float>(edge_threshold), scale); 00272 \endcode 00273 00274 <b> Required Interface:</b> 00275 00276 <ul> 00277 00278 <li> <TT>SrcIterator</TT> and <TT>DestIterator</TT> are models of ImageIterator 00279 <li> <TT>SrcAccessor</TT> and <TT>DestAccessor</TT> are models of StandardAccessor 00280 <li> <TT>SrcAccessor::value_type</TT> is a linear space 00281 <li> <TT>DiffusivityFunctor</TT> conforms to the requirements of 00282 \ref gradientBasedTransform(). Its range is between 0 and 1. 00283 <li> <TT>DiffusivityFunctor::value_type</TT> is an algebraic field 00284 00285 </ul> 00286 00287 <b> Precondition:</b> 00288 00289 <TT>scale > 0</TT> 00290 */ 00291 template <class SrcIterator, class SrcAccessor, 00292 class DestIterator, class DestAccessor, 00293 class DiffusivityFunc> 00294 void nonlinearDiffusion(SrcIterator sul, SrcIterator slr, SrcAccessor as, 00295 DestIterator dul, DestAccessor ad, 00296 DiffusivityFunc const & weight, double scale) 00297 { 00298 vigra_precondition(scale > 0.0, "nonlinearDiffusion(): scale must be > 0"); 00299 00300 double total_time = scale*scale/2.0; 00301 static const double time_step = 5.0; 00302 int number_of_steps = (int)(total_time / time_step); 00303 double rest_time = total_time - time_step * number_of_steps; 00304 00305 Size2D size(slr.x - sul.x, slr.y - sul.y); 00306 00307 typedef typename 00308 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00309 TmpType; 00310 typedef typename DiffusivityFunc::value_type WeightType; 00311 00312 BasicImage<TmpType> smooth1(size); 00313 BasicImage<TmpType> smooth2(size); 00314 00315 BasicImage<WeightType> weights(size); 00316 00317 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(), 00318 s2 = smooth2.upperLeft(); 00319 typename BasicImage<TmpType>::Accessor a = smooth1.accessor(); 00320 00321 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft(); 00322 typename BasicImage<WeightType>::Accessor wa = weights.accessor(); 00323 00324 gradientBasedTransform(sul, slr, as, wi, wa, weight); 00325 00326 internalNonlinearDiffusionAOSStep(sul, slr, as, wi, wa, s1, a, rest_time); 00327 00328 for(int i = 0; i < number_of_steps; ++i) 00329 { 00330 gradientBasedTransform(s1, s1+size, a, wi, wa, weight); 00331 00332 internalNonlinearDiffusionAOSStep(s1, s1+size, a, wi, wa, s2, a, time_step); 00333 00334 std::swap(s1, s2); 00335 } 00336 00337 copyImage(s1, s1+size, a, dul, ad); 00338 } 00339 00340 template <class SrcIterator, class SrcAccessor, 00341 class DestIterator, class DestAccessor, 00342 class DiffusivityFunc> 00343 inline 00344 void nonlinearDiffusion( 00345 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00346 pair<DestIterator, DestAccessor> dest, 00347 DiffusivityFunc const & weight, double scale) 00348 { 00349 nonlinearDiffusion(src.first, src.second, src.third, 00350 dest.first, dest.second, 00351 weight, scale); 00352 } 00353 00354 template <class SrcIterator, class SrcAccessor, 00355 class WeightIterator, class WeightAccessor, 00356 class DestIterator, class DestAccessor> 00357 void internalNonlinearDiffusionExplicitStep( 00358 SrcIterator sul, SrcIterator slr, SrcAccessor as, 00359 WeightIterator wul, WeightAccessor aw, 00360 DestIterator dul, DestAccessor ad, 00361 double time_step) 00362 { 00363 // use traits to determine SumType as to prevent possible overflow 00364 typedef typename 00365 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00366 SumType; 00367 00368 typedef typename 00369 NumericTraits<typename WeightAccessor::value_type>::RealPromote 00370 WeightType; 00371 00372 // calculate width and height of the image 00373 int w = slr.x - sul.x; 00374 int h = slr.y - sul.y; 00375 00376 int x,y; 00377 00378 static const Diff2D left(-1, 0); 00379 static const Diff2D right(1, 0); 00380 static const Diff2D top(0, -1); 00381 static const Diff2D bottom(0, 1); 00382 00383 WeightType gt, gb, gl, gr, g0; 00384 WeightType one = NumericTraits<WeightType>::one(); 00385 SumType sum; 00386 00387 time_step /= 2.0; 00388 00389 // create y iterators 00390 SrcIterator ys = sul; 00391 WeightIterator yw = wul; 00392 DestIterator yd = dul; 00393 00394 SrcIterator xs = ys; 00395 WeightIterator xw = yw; 00396 DestIterator xd = yd; 00397 00398 gt = (aw(xw) + aw(xw, bottom)) * time_step; 00399 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00400 gl = (aw(xw) + aw(xw, right)) * time_step; 00401 gr = (aw(xw) + aw(xw, right)) * time_step; 00402 g0 = one - gt - gb - gr - gl; 00403 00404 sum = g0 * as(xs); 00405 sum += gt * as(xs, bottom); 00406 sum += gb * as(xs, bottom); 00407 sum += gl * as(xs, right); 00408 sum += gr * as(xs, right); 00409 00410 ad.set(sum, xd); 00411 00412 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x) 00413 { 00414 gt = (aw(xw) + aw(xw, bottom)) * time_step; 00415 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00416 gl = (aw(xw) + aw(xw, left)) * time_step; 00417 gr = (aw(xw) + aw(xw, right)) * time_step; 00418 g0 = one - gt - gb - gr - gl; 00419 00420 sum = g0 * as(xs); 00421 sum += gt * as(xs, bottom); 00422 sum += gb * as(xs, bottom); 00423 sum += gl * as(xs, left); 00424 sum += gr * as(xs, right); 00425 00426 ad.set(sum, xd); 00427 } 00428 00429 gt = (aw(xw) + aw(xw, bottom)) * time_step; 00430 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00431 gl = (aw(xw) + aw(xw, left)) * time_step; 00432 gr = (aw(xw) + aw(xw, left)) * time_step; 00433 g0 = one - gt - gb - gr - gl; 00434 00435 sum = g0 * as(xs); 00436 sum += gt * as(xs, bottom); 00437 sum += gb * as(xs, bottom); 00438 sum += gl * as(xs, left); 00439 sum += gr * as(xs, left); 00440 00441 ad.set(sum, xd); 00442 00443 for(y=2, ++ys.y, ++yd.y, ++yw.y; y<h; ++y, ++ys.y, ++yd.y, ++yw.y) 00444 { 00445 xs = ys; 00446 xd = yd; 00447 xw = yw; 00448 00449 gt = (aw(xw) + aw(xw, top)) * time_step; 00450 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00451 gl = (aw(xw) + aw(xw, right)) * time_step; 00452 gr = (aw(xw) + aw(xw, right)) * time_step; 00453 g0 = one - gt - gb - gr - gl; 00454 00455 sum = g0 * as(xs); 00456 sum += gt * as(xs, top); 00457 sum += gb * as(xs, bottom); 00458 sum += gl * as(xs, right); 00459 sum += gr * as(xs, right); 00460 00461 ad.set(sum, xd); 00462 00463 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x) 00464 { 00465 gt = (aw(xw) + aw(xw, top)) * time_step; 00466 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00467 gl = (aw(xw) + aw(xw, left)) * time_step; 00468 gr = (aw(xw) + aw(xw, right)) * time_step; 00469 g0 = one - gt - gb - gr - gl; 00470 00471 sum = g0 * as(xs); 00472 sum += gt * as(xs, top); 00473 sum += gb * as(xs, bottom); 00474 sum += gl * as(xs, left); 00475 sum += gr * as(xs, right); 00476 00477 ad.set(sum, xd); 00478 } 00479 00480 gt = (aw(xw) + aw(xw, top)) * time_step; 00481 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00482 gl = (aw(xw) + aw(xw, left)) * time_step; 00483 gr = (aw(xw) + aw(xw, left)) * time_step; 00484 g0 = one - gt - gb - gr - gl; 00485 00486 sum = g0 * as(xs); 00487 sum += gt * as(xs, top); 00488 sum += gb * as(xs, bottom); 00489 sum += gl * as(xs, left); 00490 sum += gr * as(xs, left); 00491 00492 ad.set(sum, xd); 00493 } 00494 00495 xs = ys; 00496 xd = yd; 00497 xw = yw; 00498 00499 gt = (aw(xw) + aw(xw, top)) * time_step; 00500 gb = (aw(xw) + aw(xw, top)) * time_step; 00501 gl = (aw(xw) + aw(xw, right)) * time_step; 00502 gr = (aw(xw) + aw(xw, right)) * time_step; 00503 g0 = one - gt - gb - gr - gl; 00504 00505 sum = g0 * as(xs); 00506 sum += gt * as(xs, top); 00507 sum += gb * as(xs, top); 00508 sum += gl * as(xs, right); 00509 sum += gr * as(xs, right); 00510 00511 ad.set(sum, xd); 00512 00513 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x) 00514 { 00515 gt = (aw(xw) + aw(xw, top)) * time_step; 00516 gb = (aw(xw) + aw(xw, top)) * time_step; 00517 gl = (aw(xw) + aw(xw, left)) * time_step; 00518 gr = (aw(xw) + aw(xw, right)) * time_step; 00519 g0 = one - gt - gb - gr - gl; 00520 00521 sum = g0 * as(xs); 00522 sum += gt * as(xs, top); 00523 sum += gb * as(xs, top); 00524 sum += gl * as(xs, left); 00525 sum += gr * as(xs, right); 00526 00527 ad.set(sum, xd); 00528 } 00529 00530 gt = (aw(xw) + aw(xw, top)) * time_step; 00531 gb = (aw(xw) + aw(xw, top)) * time_step; 00532 gl = (aw(xw) + aw(xw, left)) * time_step; 00533 gr = (aw(xw) + aw(xw, left)) * time_step; 00534 g0 = one - gt - gb - gr - gl; 00535 00536 sum = g0 * as(xs); 00537 sum += gt * as(xs, top); 00538 sum += gb * as(xs, top); 00539 sum += gl * as(xs, left); 00540 sum += gr * as(xs, left); 00541 00542 ad.set(sum, xd); 00543 } 00544 00545 template <class SrcIterator, class SrcAccessor, 00546 class DestIterator, class DestAccessor, 00547 class DiffusivityFunc> 00548 void nonlinearDiffusionExplicit(SrcIterator sul, SrcIterator slr, SrcAccessor as, 00549 DestIterator dul, DestAccessor ad, 00550 DiffusivityFunc const & weight, double scale) 00551 { 00552 vigra_precondition(scale > 0.0, "nonlinearDiffusionExplicit(): scale must be > 0"); 00553 00554 double total_time = scale*scale/2.0; 00555 static const double time_step = 0.25; 00556 int number_of_steps = total_time / time_step; 00557 double rest_time = total_time - time_step * number_of_steps; 00558 00559 Size2D size(slr.x - sul.x, slr.y - sul.y); 00560 00561 typedef typename 00562 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00563 TmpType; 00564 typedef typename DiffusivityFunc::value_type WeightType; 00565 00566 BasicImage<TmpType> smooth1(size); 00567 BasicImage<TmpType> smooth2(size); 00568 00569 BasicImage<WeightType> weights(size); 00570 00571 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(), 00572 s2 = smooth2.upperLeft(); 00573 typename BasicImage<TmpType>::Accessor a = smooth1.accessor(); 00574 00575 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft(); 00576 typename BasicImage<WeightType>::Accessor wa = weights.accessor(); 00577 00578 gradientBasedTransform(sul, slr, as, wi, wa, weight); 00579 00580 internalNonlinearDiffusionExplicitStep(sul, slr, as, wi, wa, s1, a, rest_time); 00581 00582 for(int i = 0; i < number_of_steps; ++i) 00583 { 00584 gradientBasedTransform(s1, s1+size, a, wi, wa, weight); 00585 00586 internalNonlinearDiffusionExplicitStep(s1, s1+size, a, wi, wa, s2, a, time_step); 00587 00588 swap(s1, s2); 00589 } 00590 00591 copyImage(s1, s1+size, a, dul, ad); 00592 } 00593 00594 template <class SrcIterator, class SrcAccessor, 00595 class DestIterator, class DestAccessor, 00596 class DiffusivityFunc> 00597 inline 00598 void nonlinearDiffusionExplicit( 00599 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00600 pair<DestIterator, DestAccessor> dest, 00601 DiffusivityFunc const & weight, double scale) 00602 { 00603 nonlinearDiffusionExplicit(src.first, src.second, src.third, 00604 dest.first, dest.second, 00605 weight, scale); 00606 } 00607 00608 /********************************************************/ 00609 /* */ 00610 /* DiffusivityFunctor */ 00611 /* */ 00612 /********************************************************/ 00613 00614 /** \brief Diffusivity functor for non-linear diffusion. 00615 00616 This functor implements the diffusivity recommended by Weickert: 00617 00618 \f[ 00619 g(|\nabla u|) = 1 - 00620 \exp{\left(\frac{-3.315}{(|\nabla u| / thresh)^4}\right)} 00621 \f] 00622 00623 00624 where <TT>thresh</TT> is a threshold that determines whether a specific gradient 00625 magnitude is interpreted as a significant edge (above the threshold) 00626 or as noise. It is meant to be used with \ref nonlinearDiffusion(). 00627 */ 00628 template <class Value> 00629 class DiffusivityFunctor 00630 { 00631 public: 00632 /** the functors first argument type (must be an algebraic field with <TT>exp()</TT> defined). 00633 However, the functor also works with RGBValue<first_argument_type> (this hack is 00634 necessary since Microsoft C++ doesn't support partial specialization). 00635 */ 00636 typedef Value first_argument_type; 00637 00638 /** the functors second argument type (same as first). 00639 However, the functor also works with RGBValue<second_argument_type> (this hack is 00640 necessary since Microsoft C++ doesn't support partial specialization). 00641 */ 00642 typedef Value second_argument_type; 00643 00644 /** the functors result type 00645 */ 00646 typedef typename NumericTraits<Value>::RealPromote result_type; 00647 00648 /** \deprecated use first_argument_type, second_argument_type, result_type 00649 */ 00650 typedef Value value_type; 00651 00652 /** init functor with given edge threshold 00653 */ 00654 DiffusivityFunctor(Value const & thresh) 00655 : weight_(thresh*thresh), 00656 one_(NumericTraits<result_type>::one()), 00657 zero_(NumericTraits<result_type>::zero()) 00658 {} 00659 00660 /** calculate diffusivity from scalar arguments 00661 */ 00662 result_type operator()(first_argument_type const & gx, second_argument_type const & gy) const 00663 { 00664 Value mag = (gx*gx + gy*gy) / weight_; 00665 00666 return (mag == zero_) ? one_ : one_ - exp(-3.315 / mag / mag); 00667 } 00668 00669 /** calculate diffusivity from RGB arguments 00670 */ 00671 result_type operator()(RGBValue<Value> const & gx, RGBValue<Value> const & gy) const 00672 { 00673 result_type mag = (gx.red()*gx.red() + 00674 gx.green()*gx.green() + 00675 gx.blue()*gx.blue() + 00676 gy.red()*gy.red() + 00677 gy.green()*gy.green() + 00678 gy.blue()*gy.blue()) / weight_; 00679 00680 return (mag == zero_) ? one_ : one_ - exp(-3.315 / mag / mag); 00681 } 00682 00683 result_type weight_; 00684 result_type one_; 00685 result_type zero_; 00686 }; 00687 00688 //@} 00689 00690 } // namespace vigra 00691 00692 #endif /* VIGRA_NONLINEARDIFFUSION_HXX */
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|