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