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