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