[ 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 /* ( Version 1.6.0, Aug 13 2008 ) */ 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 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00012 /* vigra@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 doxygen_overloaded_function(template <...> void nonlinearDiffusion) 00308 00309 template <class SrcIterator, class SrcAccessor, 00310 class DestIterator, class DestAccessor, 00311 class DiffusivityFunc> 00312 void nonlinearDiffusion(SrcIterator sul, SrcIterator slr, SrcAccessor as, 00313 DestIterator dul, DestAccessor ad, 00314 DiffusivityFunc const & weight, double scale) 00315 { 00316 vigra_precondition(scale > 0.0, "nonlinearDiffusion(): scale must be > 0"); 00317 00318 double total_time = scale*scale/2.0; 00319 static const double time_step = 5.0; 00320 int number_of_steps = (int)(total_time / time_step); 00321 double rest_time = total_time - time_step * number_of_steps; 00322 00323 Size2D size(slr.x - sul.x, slr.y - sul.y); 00324 00325 typedef typename 00326 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00327 TmpType; 00328 typedef typename DiffusivityFunc::value_type WeightType; 00329 00330 BasicImage<TmpType> smooth1(size); 00331 BasicImage<TmpType> smooth2(size); 00332 00333 BasicImage<WeightType> weights(size); 00334 00335 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(), 00336 s2 = smooth2.upperLeft(); 00337 typename BasicImage<TmpType>::Accessor a = smooth1.accessor(); 00338 00339 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft(); 00340 typename BasicImage<WeightType>::Accessor wa = weights.accessor(); 00341 00342 gradientBasedTransform(sul, slr, as, wi, wa, weight); 00343 00344 internalNonlinearDiffusionAOSStep(sul, slr, as, wi, wa, s1, a, rest_time); 00345 00346 for(int i = 0; i < number_of_steps; ++i) 00347 { 00348 gradientBasedTransform(s1, s1+size, a, wi, wa, weight); 00349 00350 internalNonlinearDiffusionAOSStep(s1, s1+size, a, wi, wa, s2, a, time_step); 00351 00352 std::swap(s1, s2); 00353 } 00354 00355 copyImage(s1, s1+size, a, dul, ad); 00356 } 00357 00358 template <class SrcIterator, class SrcAccessor, 00359 class DestIterator, class DestAccessor, 00360 class DiffusivityFunc> 00361 inline 00362 void nonlinearDiffusion( 00363 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00364 pair<DestIterator, DestAccessor> dest, 00365 DiffusivityFunc const & weight, double scale) 00366 { 00367 nonlinearDiffusion(src.first, src.second, src.third, 00368 dest.first, dest.second, 00369 weight, scale); 00370 } 00371 00372 template <class SrcIterator, class SrcAccessor, 00373 class WeightIterator, class WeightAccessor, 00374 class DestIterator, class DestAccessor> 00375 void internalNonlinearDiffusionExplicitStep( 00376 SrcIterator sul, SrcIterator slr, SrcAccessor as, 00377 WeightIterator wul, WeightAccessor aw, 00378 DestIterator dul, DestAccessor ad, 00379 double time_step) 00380 { 00381 // use traits to determine SumType as to prevent possible overflow 00382 typedef typename 00383 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00384 SumType; 00385 00386 typedef typename 00387 NumericTraits<typename WeightAccessor::value_type>::RealPromote 00388 WeightType; 00389 00390 // calculate width and height of the image 00391 int w = slr.x - sul.x; 00392 int h = slr.y - sul.y; 00393 00394 int x,y; 00395 00396 static const Diff2D left(-1, 0); 00397 static const Diff2D right(1, 0); 00398 static const Diff2D top(0, -1); 00399 static const Diff2D bottom(0, 1); 00400 00401 WeightType gt, gb, gl, gr, g0; 00402 WeightType one = NumericTraits<WeightType>::one(); 00403 SumType sum; 00404 00405 time_step /= 2.0; 00406 00407 // create y iterators 00408 SrcIterator ys = sul; 00409 WeightIterator yw = wul; 00410 DestIterator yd = dul; 00411 00412 SrcIterator xs = ys; 00413 WeightIterator xw = yw; 00414 DestIterator xd = yd; 00415 00416 gt = (aw(xw) + aw(xw, bottom)) * time_step; 00417 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00418 gl = (aw(xw) + aw(xw, right)) * time_step; 00419 gr = (aw(xw) + aw(xw, right)) * time_step; 00420 g0 = one - gt - gb - gr - gl; 00421 00422 sum = g0 * as(xs); 00423 sum += gt * as(xs, bottom); 00424 sum += gb * as(xs, bottom); 00425 sum += gl * as(xs, right); 00426 sum += gr * as(xs, right); 00427 00428 ad.set(sum, xd); 00429 00430 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x) 00431 { 00432 gt = (aw(xw) + aw(xw, bottom)) * time_step; 00433 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00434 gl = (aw(xw) + aw(xw, left)) * time_step; 00435 gr = (aw(xw) + aw(xw, right)) * time_step; 00436 g0 = one - gt - gb - gr - gl; 00437 00438 sum = g0 * as(xs); 00439 sum += gt * as(xs, bottom); 00440 sum += gb * as(xs, bottom); 00441 sum += gl * as(xs, left); 00442 sum += gr * as(xs, right); 00443 00444 ad.set(sum, xd); 00445 } 00446 00447 gt = (aw(xw) + aw(xw, bottom)) * time_step; 00448 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00449 gl = (aw(xw) + aw(xw, left)) * time_step; 00450 gr = (aw(xw) + aw(xw, left)) * time_step; 00451 g0 = one - gt - gb - gr - gl; 00452 00453 sum = g0 * as(xs); 00454 sum += gt * as(xs, bottom); 00455 sum += gb * as(xs, bottom); 00456 sum += gl * as(xs, left); 00457 sum += gr * as(xs, left); 00458 00459 ad.set(sum, xd); 00460 00461 for(y=2, ++ys.y, ++yd.y, ++yw.y; y<h; ++y, ++ys.y, ++yd.y, ++yw.y) 00462 { 00463 xs = ys; 00464 xd = yd; 00465 xw = yw; 00466 00467 gt = (aw(xw) + aw(xw, top)) * time_step; 00468 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00469 gl = (aw(xw) + aw(xw, right)) * time_step; 00470 gr = (aw(xw) + aw(xw, right)) * time_step; 00471 g0 = one - gt - gb - gr - gl; 00472 00473 sum = g0 * as(xs); 00474 sum += gt * as(xs, top); 00475 sum += gb * as(xs, bottom); 00476 sum += gl * as(xs, right); 00477 sum += gr * as(xs, right); 00478 00479 ad.set(sum, xd); 00480 00481 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x) 00482 { 00483 gt = (aw(xw) + aw(xw, top)) * time_step; 00484 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00485 gl = (aw(xw) + aw(xw, left)) * time_step; 00486 gr = (aw(xw) + aw(xw, right)) * time_step; 00487 g0 = one - gt - gb - gr - gl; 00488 00489 sum = g0 * as(xs); 00490 sum += gt * as(xs, top); 00491 sum += gb * as(xs, bottom); 00492 sum += gl * as(xs, left); 00493 sum += gr * as(xs, right); 00494 00495 ad.set(sum, xd); 00496 } 00497 00498 gt = (aw(xw) + aw(xw, top)) * time_step; 00499 gb = (aw(xw) + aw(xw, bottom)) * time_step; 00500 gl = (aw(xw) + aw(xw, left)) * time_step; 00501 gr = (aw(xw) + aw(xw, left)) * time_step; 00502 g0 = one - gt - gb - gr - gl; 00503 00504 sum = g0 * as(xs); 00505 sum += gt * as(xs, top); 00506 sum += gb * as(xs, bottom); 00507 sum += gl * as(xs, left); 00508 sum += gr * as(xs, left); 00509 00510 ad.set(sum, xd); 00511 } 00512 00513 xs = ys; 00514 xd = yd; 00515 xw = yw; 00516 00517 gt = (aw(xw) + aw(xw, top)) * time_step; 00518 gb = (aw(xw) + aw(xw, top)) * time_step; 00519 gl = (aw(xw) + aw(xw, right)) * time_step; 00520 gr = (aw(xw) + aw(xw, right)) * time_step; 00521 g0 = one - gt - gb - gr - gl; 00522 00523 sum = g0 * as(xs); 00524 sum += gt * as(xs, top); 00525 sum += gb * as(xs, top); 00526 sum += gl * as(xs, right); 00527 sum += gr * as(xs, right); 00528 00529 ad.set(sum, xd); 00530 00531 for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x) 00532 { 00533 gt = (aw(xw) + aw(xw, top)) * time_step; 00534 gb = (aw(xw) + aw(xw, top)) * time_step; 00535 gl = (aw(xw) + aw(xw, left)) * time_step; 00536 gr = (aw(xw) + aw(xw, right)) * time_step; 00537 g0 = one - gt - gb - gr - gl; 00538 00539 sum = g0 * as(xs); 00540 sum += gt * as(xs, top); 00541 sum += gb * as(xs, top); 00542 sum += gl * as(xs, left); 00543 sum += gr * as(xs, right); 00544 00545 ad.set(sum, xd); 00546 } 00547 00548 gt = (aw(xw) + aw(xw, top)) * time_step; 00549 gb = (aw(xw) + aw(xw, top)) * time_step; 00550 gl = (aw(xw) + aw(xw, left)) * time_step; 00551 gr = (aw(xw) + aw(xw, left)) * time_step; 00552 g0 = one - gt - gb - gr - gl; 00553 00554 sum = g0 * as(xs); 00555 sum += gt * as(xs, top); 00556 sum += gb * as(xs, top); 00557 sum += gl * as(xs, left); 00558 sum += gr * as(xs, left); 00559 00560 ad.set(sum, xd); 00561 } 00562 00563 template <class SrcIterator, class SrcAccessor, 00564 class DestIterator, class DestAccessor, 00565 class DiffusivityFunc> 00566 void nonlinearDiffusionExplicit(SrcIterator sul, SrcIterator slr, SrcAccessor as, 00567 DestIterator dul, DestAccessor ad, 00568 DiffusivityFunc const & weight, double scale) 00569 { 00570 vigra_precondition(scale > 0.0, "nonlinearDiffusionExplicit(): scale must be > 0"); 00571 00572 double total_time = scale*scale/2.0; 00573 static const double time_step = 0.25; 00574 int number_of_steps = total_time / time_step; 00575 double rest_time = total_time - time_step * number_of_steps; 00576 00577 Size2D size(slr.x - sul.x, slr.y - sul.y); 00578 00579 typedef typename 00580 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00581 TmpType; 00582 typedef typename DiffusivityFunc::value_type WeightType; 00583 00584 BasicImage<TmpType> smooth1(size); 00585 BasicImage<TmpType> smooth2(size); 00586 00587 BasicImage<WeightType> weights(size); 00588 00589 typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(), 00590 s2 = smooth2.upperLeft(); 00591 typename BasicImage<TmpType>::Accessor a = smooth1.accessor(); 00592 00593 typename BasicImage<WeightType>::Iterator wi = weights.upperLeft(); 00594 typename BasicImage<WeightType>::Accessor wa = weights.accessor(); 00595 00596 gradientBasedTransform(sul, slr, as, wi, wa, weight); 00597 00598 internalNonlinearDiffusionExplicitStep(sul, slr, as, wi, wa, s1, a, rest_time); 00599 00600 for(int i = 0; i < number_of_steps; ++i) 00601 { 00602 gradientBasedTransform(s1, s1+size, a, wi, wa, weight); 00603 00604 internalNonlinearDiffusionExplicitStep(s1, s1+size, a, wi, wa, s2, a, time_step); 00605 00606 swap(s1, s2); 00607 } 00608 00609 copyImage(s1, s1+size, a, dul, ad); 00610 } 00611 00612 template <class SrcIterator, class SrcAccessor, 00613 class DestIterator, class DestAccessor, 00614 class DiffusivityFunc> 00615 inline 00616 void nonlinearDiffusionExplicit( 00617 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00618 pair<DestIterator, DestAccessor> dest, 00619 DiffusivityFunc const & weight, double scale) 00620 { 00621 nonlinearDiffusionExplicit(src.first, src.second, src.third, 00622 dest.first, dest.second, 00623 weight, scale); 00624 } 00625 00626 /********************************************************/ 00627 /* */ 00628 /* DiffusivityFunctor */ 00629 /* */ 00630 /********************************************************/ 00631 00632 /** \brief Diffusivity functor for non-linear diffusion. 00633 00634 This functor implements the diffusivity recommended by Weickert: 00635 00636 \f[ 00637 g(|\nabla u|) = 1 - 00638 \exp{\left(\frac{-3.315}{(|\nabla u| / thresh)^4}\right)} 00639 \f] 00640 00641 00642 where <TT>thresh</TT> is a threshold that determines whether a specific gradient 00643 magnitude is interpreted as a significant edge (above the threshold) 00644 or as noise. It is meant to be used with \ref nonlinearDiffusion(). 00645 */ 00646 template <class Value> 00647 class DiffusivityFunctor 00648 { 00649 public: 00650 /** the functors first argument type (must be an algebraic field with <TT>exp()</TT> defined). 00651 However, the functor also works with RGBValue<first_argument_type> (this hack is 00652 necessary since Microsoft C++ doesn't support partial specialization). 00653 */ 00654 typedef Value first_argument_type; 00655 00656 /** the functors second argument type (same as first). 00657 However, the functor also works with RGBValue<second_argument_type> (this hack is 00658 necessary since Microsoft C++ doesn't support partial specialization). 00659 */ 00660 typedef Value second_argument_type; 00661 00662 /** the functors result type 00663 */ 00664 typedef typename NumericTraits<Value>::RealPromote result_type; 00665 00666 /** \deprecated use first_argument_type, second_argument_type, result_type 00667 */ 00668 typedef Value value_type; 00669 00670 /** init functor with given edge threshold 00671 */ 00672 DiffusivityFunctor(Value const & thresh) 00673 : weight_(thresh*thresh), 00674 one_(NumericTraits<result_type>::one()), 00675 zero_(NumericTraits<result_type>::zero()) 00676 {} 00677 00678 /** calculate diffusivity from scalar arguments 00679 */ 00680 result_type operator()(first_argument_type const & gx, second_argument_type const & gy) const 00681 { 00682 Value mag = (gx*gx + gy*gy) / weight_; 00683 00684 return (mag == zero_) ? one_ : one_ - VIGRA_CSTD::exp(-3.315 / mag / mag); 00685 } 00686 00687 /** calculate diffusivity from RGB arguments 00688 */ 00689 result_type operator()(RGBValue<Value> const & gx, RGBValue<Value> const & gy) const 00690 { 00691 result_type mag = (gx.red()*gx.red() + 00692 gx.green()*gx.green() + 00693 gx.blue()*gx.blue() + 00694 gy.red()*gy.red() + 00695 gy.green()*gy.green() + 00696 gy.blue()*gy.blue()) / weight_; 00697 00698 return (mag == zero_) ? one_ : one_ - VIGRA_CSTD::exp(-3.315 / mag / mag); 00699 } 00700 00701 result_type weight_; 00702 result_type one_; 00703 result_type zero_; 00704 }; 00705 00706 template <class ValueType> 00707 class FunctorTraits<DiffusivityFunctor<ValueType> > 00708 : public FunctorTraitsBase<DiffusivityFunctor<ValueType> > 00709 { 00710 public: 00711 typedef VigraTrueType isBinaryFunctor; 00712 }; 00713 00714 00715 //@} 00716 00717 } // namespace vigra 00718 00719 #endif /* VIGRA_NONLINEARDIFFUSION_HXX */
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|