vigra - vigra/nonlineardiffusion.hxx Source File

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

vigra/nonlineardiffusion.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)