[ 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 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@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 "stdimage.hxx"
00043 #include "stdimagefunctions.hxx"
00044 #include "imageiteratoradapter.hxx"
00045 #include "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 doxygen_overloaded_function(template <...> void nonlinearDiffusion)
00308 
00309 template <class SrcIterator, class SrcAccessor,
00310           class DestIterator, class DestAccessor,
00311           class DiffusivityFunc>
00312 void nonlinearDiffusion(SrcIterator sul, SrcIterator slr, SrcAccessor as,
00313                    DestIterator dul, DestAccessor ad,
00314                    DiffusivityFunc const & weight, double scale)
00315 {
00316     vigra_precondition(scale > 0.0, "nonlinearDiffusion(): scale must be > 0");
00317     
00318     double total_time = scale*scale/2.0;
00319     static const double time_step = 5.0;
00320     int number_of_steps = (int)(total_time / time_step);
00321     double rest_time = total_time - time_step * number_of_steps;
00322     
00323     Size2D size(slr.x - sul.x, slr.y - sul.y);
00324 
00325     typedef typename
00326         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00327         TmpType;
00328     typedef typename DiffusivityFunc::value_type WeightType;
00329     
00330     BasicImage<TmpType> smooth1(size);
00331     BasicImage<TmpType> smooth2(size);
00332     
00333     BasicImage<WeightType> weights(size);
00334     
00335     typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
00336                                   s2 = smooth2.upperLeft();
00337     typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
00338     
00339     typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
00340     typename BasicImage<WeightType>::Accessor wa = weights.accessor();
00341 
00342     gradientBasedTransform(sul, slr, as, wi, wa, weight);
00343 
00344     internalNonlinearDiffusionAOSStep(sul, slr, as, wi, wa, s1, a, rest_time);
00345 
00346     for(int i = 0; i < number_of_steps; ++i)
00347     {
00348         gradientBasedTransform(s1, s1+size, a, wi, wa, weight);
00349                       
00350         internalNonlinearDiffusionAOSStep(s1, s1+size, a, wi, wa, s2, a, time_step);
00351     
00352         std::swap(s1, s2);
00353     }
00354     
00355     copyImage(s1, s1+size, a, dul, ad);
00356 }
00357 
00358 template <class SrcIterator, class SrcAccessor,
00359           class DestIterator, class DestAccessor,
00360           class DiffusivityFunc>
00361 inline
00362 void nonlinearDiffusion(
00363     triple<SrcIterator, SrcIterator, SrcAccessor> src,
00364     pair<DestIterator, DestAccessor> dest,
00365     DiffusivityFunc const & weight, double scale)
00366 {
00367     nonlinearDiffusion(src.first, src.second, src.third,
00368                            dest.first, dest.second,
00369                            weight, scale);
00370 }
00371 
00372 template <class SrcIterator, class SrcAccessor,
00373           class WeightIterator, class WeightAccessor,
00374           class DestIterator, class DestAccessor>
00375 void internalNonlinearDiffusionExplicitStep(
00376                    SrcIterator sul, SrcIterator slr, SrcAccessor as,
00377                    WeightIterator wul, WeightAccessor aw,
00378                    DestIterator dul, DestAccessor ad,
00379                    double time_step)
00380 {
00381     // use traits to determine SumType as to prevent possible overflow
00382     typedef typename
00383         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00384         SumType;
00385     
00386     typedef typename
00387         NumericTraits<typename WeightAccessor::value_type>::RealPromote
00388         WeightType;
00389         
00390     // calculate width and height of the image
00391     int w = slr.x - sul.x;
00392     int h = slr.y - sul.y;
00393 
00394     int x,y;
00395     
00396     static const Diff2D left(-1, 0);
00397     static const Diff2D right(1, 0);
00398     static const Diff2D top(0, -1);
00399     static const Diff2D bottom(0, 1);
00400     
00401     WeightType gt, gb, gl, gr, g0;
00402     WeightType one = NumericTraits<WeightType>::one();
00403     SumType sum;
00404     
00405     time_step /= 2.0;
00406     
00407     // create y iterators
00408     SrcIterator ys = sul;
00409     WeightIterator yw = wul;
00410     DestIterator yd = dul;
00411         
00412     SrcIterator xs = ys;
00413     WeightIterator xw = yw;
00414     DestIterator xd = yd;
00415     
00416     gt = (aw(xw) + aw(xw, bottom)) * time_step;
00417     gb = (aw(xw) + aw(xw, bottom)) * time_step;
00418     gl = (aw(xw) + aw(xw, right)) * time_step;
00419     gr = (aw(xw) + aw(xw, right)) * time_step;
00420     g0 = one - gt - gb - gr - gl;
00421 
00422     sum = g0 * as(xs);
00423     sum += gt * as(xs, bottom);
00424     sum += gb * as(xs, bottom);
00425     sum += gl * as(xs, right);
00426     sum += gr * as(xs, right);
00427 
00428     ad.set(sum, xd);
00429 
00430     for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
00431     {
00432         gt = (aw(xw) + aw(xw, bottom)) * time_step;
00433         gb = (aw(xw) + aw(xw, bottom)) * time_step;
00434         gl = (aw(xw) + aw(xw, left)) * time_step;
00435         gr = (aw(xw) + aw(xw, right)) * time_step;
00436         g0 = one - gt - gb - gr - gl;
00437 
00438         sum = g0 * as(xs);
00439         sum += gt * as(xs, bottom);
00440         sum += gb * as(xs, bottom);
00441         sum += gl * as(xs, left);
00442         sum += gr * as(xs, right);
00443 
00444         ad.set(sum, xd);
00445     }
00446 
00447     gt = (aw(xw) + aw(xw, bottom)) * time_step;
00448     gb = (aw(xw) + aw(xw, bottom)) * time_step;
00449     gl = (aw(xw) + aw(xw, left)) * time_step;
00450     gr = (aw(xw) + aw(xw, left)) * time_step;
00451     g0 = one - gt - gb - gr - gl;
00452 
00453     sum = g0 * as(xs);
00454     sum += gt * as(xs, bottom);
00455     sum += gb * as(xs, bottom);
00456     sum += gl * as(xs, left);
00457     sum += gr * as(xs, left);
00458 
00459     ad.set(sum, xd);
00460     
00461     for(y=2, ++ys.y, ++yd.y, ++yw.y; y<h; ++y, ++ys.y, ++yd.y, ++yw.y)
00462     {
00463         xs = ys;
00464         xd = yd;
00465         xw = yw;
00466         
00467         gt = (aw(xw) + aw(xw, top)) * time_step;
00468         gb = (aw(xw) + aw(xw, bottom)) * time_step;
00469         gl = (aw(xw) + aw(xw, right)) * time_step;
00470         gr = (aw(xw) + aw(xw, right)) * time_step;
00471         g0 = one - gt - gb - gr - gl;
00472 
00473         sum = g0 * as(xs);
00474         sum += gt * as(xs, top);
00475         sum += gb * as(xs, bottom);
00476         sum += gl * as(xs, right);
00477         sum += gr * as(xs, right);
00478 
00479         ad.set(sum, xd);
00480         
00481         for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
00482         {
00483             gt = (aw(xw) + aw(xw, top)) * time_step;
00484             gb = (aw(xw) + aw(xw, bottom)) * time_step;
00485             gl = (aw(xw) + aw(xw, left)) * time_step;
00486             gr = (aw(xw) + aw(xw, right)) * time_step;
00487             g0 = one - gt - gb - gr - gl;
00488             
00489             sum = g0 * as(xs);
00490             sum += gt * as(xs, top);
00491             sum += gb * as(xs, bottom);
00492             sum += gl * as(xs, left);
00493             sum += gr * as(xs, right);
00494             
00495             ad.set(sum, xd);
00496         }
00497         
00498         gt = (aw(xw) + aw(xw, top)) * time_step;
00499         gb = (aw(xw) + aw(xw, bottom)) * time_step;
00500         gl = (aw(xw) + aw(xw, left)) * time_step;
00501         gr = (aw(xw) + aw(xw, left)) * time_step;
00502         g0 = one - gt - gb - gr - gl;
00503 
00504         sum = g0 * as(xs);
00505         sum += gt * as(xs, top);
00506         sum += gb * as(xs, bottom);
00507         sum += gl * as(xs, left);
00508         sum += gr * as(xs, left);
00509 
00510         ad.set(sum, xd);
00511     }
00512     
00513     xs = ys;
00514     xd = yd;
00515     xw = yw;
00516 
00517     gt = (aw(xw) + aw(xw, top)) * time_step;
00518     gb = (aw(xw) + aw(xw, top)) * time_step;
00519     gl = (aw(xw) + aw(xw, right)) * time_step;
00520     gr = (aw(xw) + aw(xw, right)) * time_step;
00521     g0 = one - gt - gb - gr - gl;
00522 
00523     sum = g0 * as(xs);
00524     sum += gt * as(xs, top);
00525     sum += gb * as(xs, top);
00526     sum += gl * as(xs, right);
00527     sum += gr * as(xs, right);
00528 
00529     ad.set(sum, xd);
00530 
00531     for(x=2, ++xs.x, ++xd.x, ++xw.x; x<w; ++x, ++xs.x, ++xd.x, ++xw.x)
00532     {
00533         gt = (aw(xw) + aw(xw, top)) * time_step;
00534         gb = (aw(xw) + aw(xw, top)) * time_step;
00535         gl = (aw(xw) + aw(xw, left)) * time_step;
00536         gr = (aw(xw) + aw(xw, right)) * time_step;
00537         g0 = one - gt - gb - gr - gl;
00538 
00539         sum = g0 * as(xs);
00540         sum += gt * as(xs, top);
00541         sum += gb * as(xs, top);
00542         sum += gl * as(xs, left);
00543         sum += gr * as(xs, right);
00544 
00545         ad.set(sum, xd);
00546     }
00547 
00548     gt = (aw(xw) + aw(xw, top)) * time_step;
00549     gb = (aw(xw) + aw(xw, top)) * time_step;
00550     gl = (aw(xw) + aw(xw, left)) * time_step;
00551     gr = (aw(xw) + aw(xw, left)) * time_step;
00552     g0 = one - gt - gb - gr - gl;
00553 
00554     sum = g0 * as(xs);
00555     sum += gt * as(xs, top);
00556     sum += gb * as(xs, top);
00557     sum += gl * as(xs, left);
00558     sum += gr * as(xs, left);
00559 
00560     ad.set(sum, xd);
00561 }
00562 
00563 template <class SrcIterator, class SrcAccessor,
00564           class DestIterator, class DestAccessor,
00565           class DiffusivityFunc>
00566 void nonlinearDiffusionExplicit(SrcIterator sul, SrcIterator slr, SrcAccessor as,
00567                    DestIterator dul, DestAccessor ad,
00568                    DiffusivityFunc const & weight, double scale)
00569 {
00570     vigra_precondition(scale > 0.0, "nonlinearDiffusionExplicit(): scale must be > 0");
00571     
00572     double total_time = scale*scale/2.0;
00573     static const double time_step = 0.25;
00574     int number_of_steps = total_time / time_step;
00575     double rest_time = total_time - time_step * number_of_steps;
00576     
00577     Size2D size(slr.x - sul.x, slr.y - sul.y);
00578 
00579     typedef typename
00580         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00581         TmpType;
00582     typedef typename DiffusivityFunc::value_type WeightType;
00583     
00584     BasicImage<TmpType> smooth1(size);
00585     BasicImage<TmpType> smooth2(size);
00586     
00587     BasicImage<WeightType> weights(size);
00588     
00589     typename BasicImage<TmpType>::Iterator s1 = smooth1.upperLeft(),
00590                                   s2 = smooth2.upperLeft();
00591     typename BasicImage<TmpType>::Accessor a = smooth1.accessor();
00592     
00593     typename BasicImage<WeightType>::Iterator wi = weights.upperLeft();
00594     typename BasicImage<WeightType>::Accessor wa = weights.accessor();
00595 
00596     gradientBasedTransform(sul, slr, as, wi, wa, weight);
00597 
00598     internalNonlinearDiffusionExplicitStep(sul, slr, as, wi, wa, s1, a, rest_time);
00599 
00600     for(int i = 0; i < number_of_steps; ++i)
00601     {
00602         gradientBasedTransform(s1, s1+size, a, wi, wa, weight);
00603                       
00604         internalNonlinearDiffusionExplicitStep(s1, s1+size, a, wi, wa, s2, a, time_step);
00605     
00606         swap(s1, s2);
00607     }
00608     
00609     copyImage(s1, s1+size, a, dul, ad);
00610 }
00611 
00612 template <class SrcIterator, class SrcAccessor,
00613           class DestIterator, class DestAccessor,
00614           class DiffusivityFunc>
00615 inline
00616 void nonlinearDiffusionExplicit(
00617     triple<SrcIterator, SrcIterator, SrcAccessor> src,
00618     pair<DestIterator, DestAccessor> dest,
00619     DiffusivityFunc const & weight, double scale)
00620 {
00621     nonlinearDiffusionExplicit(src.first, src.second, src.third,
00622                            dest.first, dest.second,
00623                            weight, scale);
00624 }
00625 
00626 /********************************************************/
00627 /*                                                      */
00628 /*                   DiffusivityFunctor                 */
00629 /*                                                      */
00630 /********************************************************/
00631 
00632 /** \brief Diffusivity functor for non-linear diffusion.
00633 
00634     This functor implements the diffusivity recommended by Weickert:
00635     
00636     \f[
00637         g(|\nabla u|) = 1 -
00638            \exp{\left(\frac{-3.315}{(|\nabla u| / thresh)^4}\right)}
00639     \f]
00640     
00641     
00642     where <TT>thresh</TT> is a threshold that determines whether a specific gradient
00643     magnitude is interpreted as a significant edge (above the threshold)
00644     or as noise. It is meant to be used with \ref nonlinearDiffusion().
00645 */
00646 template <class Value>
00647 class DiffusivityFunctor
00648 {
00649   public:
00650          /** the functors first argument type (must be an algebraic field with <TT>exp()</TT> defined).
00651              However, the functor also works with RGBValue<first_argument_type> (this hack is
00652              necessary since Microsoft C++ doesn't support partial specialization).
00653          */
00654     typedef Value first_argument_type;
00655     
00656          /** the functors second argument type (same as first).
00657              However, the functor also works with RGBValue<second_argument_type> (this hack is
00658              necessary since Microsoft C++ doesn't support partial specialization).
00659          */
00660     typedef Value second_argument_type;
00661     
00662          /** the functors result type
00663          */
00664     typedef typename NumericTraits<Value>::RealPromote result_type;
00665     
00666          /** \deprecated use first_argument_type, second_argument_type, result_type
00667          */
00668     typedef Value value_type;
00669     
00670          /** init functor with given edge threshold
00671          */
00672     DiffusivityFunctor(Value const & thresh)
00673     : weight_(thresh*thresh),
00674       one_(NumericTraits<result_type>::one()),
00675       zero_(NumericTraits<result_type>::zero())
00676     {}
00677     
00678          /** calculate diffusivity from scalar arguments
00679          */
00680     result_type operator()(first_argument_type const & gx, second_argument_type const & gy) const
00681     {
00682         Value mag = (gx*gx + gy*gy) / weight_;
00683                      
00684         return (mag == zero_) ? one_ : one_ - VIGRA_CSTD::exp(-3.315 / mag / mag);
00685     }
00686     
00687          /** calculate diffusivity from RGB arguments
00688          */
00689     result_type operator()(RGBValue<Value> const & gx, RGBValue<Value> const & gy) const
00690     {
00691         result_type mag = (gx.red()*gx.red() +
00692                      gx.green()*gx.green() +
00693                      gx.blue()*gx.blue() +
00694                      gy.red()*gy.red() +
00695                      gy.green()*gy.green() +
00696                      gy.blue()*gy.blue()) / weight_;
00697 
00698         return (mag == zero_) ? one_ : one_ - VIGRA_CSTD::exp(-3.315 / mag / mag);
00699     }
00700     
00701     result_type weight_;
00702     result_type one_;
00703     result_type zero_;
00704 };
00705 
00706 template <class ValueType>
00707 class FunctorTraits<DiffusivityFunctor<ValueType> >
00708 : public FunctorTraitsBase<DiffusivityFunctor<ValueType> >
00709 {
00710   public:
00711     typedef VigraTrueType isBinaryFunctor;
00712 };
00713 
00714 
00715 //@}
00716 
00717 } // namespace vigra
00718 
00719 #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 (13 Aug 2008)