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

details vigra/nonlineardiffusion.hxx VIGRA

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