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

| vigra/nonlineardiffusion.hxx |
--- | --- | ---

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.3.2, Jan 27 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.2 (27 Jan 2005)