[ 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
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     propotional 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> <<a href="nonlineardiffusion_8hxx-source.html">vigra/nonlineardiffusion.hxx</a>>
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.7.0 (Thu Aug 25 2011)