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