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