[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]

details vigra/nonlineardiffusion.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.3 (18 Aug 2005)