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