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

vigra/watersheds.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2005 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_WATERSHEDS_HXX
00037 #define VIGRA_WATERSHEDS_HXX
00038 
00039 #include <functional>
00040 #include "mathutil.hxx"
00041 #include "stdimage.hxx"
00042 #include "pixelneighborhood.hxx"
00043 #include "localminmax.hxx"
00044 #include "labelimage.hxx"
00045 #include "seededregiongrowing.hxx"
00046 #include "functorexpression.hxx"
00047 #include "union_find.hxx"
00048 
00049 namespace vigra {
00050 
00051 template <class SrcIterator, class SrcAccessor,
00052           class DestIterator, class DestAccessor,
00053           class Neighborhood>
00054 unsigned int watershedLabeling(SrcIterator upperlefts,
00055                         SrcIterator lowerrights, SrcAccessor sa,
00056                         DestIterator upperleftd, DestAccessor da,
00057                         Neighborhood)
00058 {
00059     typedef typename DestAccessor::value_type LabelType;
00060 
00061     int w = lowerrights.x - upperlefts.x;
00062     int h = lowerrights.y - upperlefts.y;
00063     int x,y;
00064 
00065     SrcIterator ys(upperlefts);
00066     SrcIterator xs(ys);
00067     DestIterator yd(upperleftd);
00068     DestIterator xd(yd);
00069 
00070     // temporary image to store region labels
00071     detail::UnionFindArray<LabelType> labels;
00072 
00073     // initialize the neighborhood circulators
00074     NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst);
00075     NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North);
00076     NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast);
00077     ++ncend;
00078     NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North);
00079     ++ncendBorder;
00080 
00081     // pass 1: scan image from upper left to lower right
00082     // to find connected components
00083 
00084     // Each component will be represented by a tree of pixels. Each
00085     // pixel contains the scan order address of its parent in the
00086     // tree.  In order for pass 2 to work correctly, the parent must
00087     // always have a smaller scan order address than the child.
00088     // Therefore, we can merge trees only at their roots, because the
00089     // root of the combined tree must have the smallest scan order
00090     // address among all the tree's pixels/ nodes.  The root of each
00091     // tree is distinguished by pointing to itself (it contains its
00092     // own scan order address). This condition is enforced whenever a
00093     // new region is found or two regions are merged
00094     da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
00095 
00096     ++xs.x;
00097     ++xd.x;
00098     for(x = 1; x != w; ++x, ++xs.x, ++xd.x)
00099     {
00100         if((sa(xs) & Neighborhood::directionBit(Neighborhood::West)) ||
00101            (sa(xs, Neighborhood::west()) & Neighborhood::directionBit(Neighborhood::East)))
00102         {
00103             da.set(da(xd, Neighborhood::west()), xd);
00104         }
00105         else
00106         {
00107             da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
00108         }
00109     }
00110 
00111     ++ys.y;
00112     ++yd.y;
00113     for(y = 1; y != h; ++y, ++ys.y, ++yd.y)
00114     {
00115         xs = ys;
00116         xd = yd;
00117 
00118         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00119         {
00120             NeighborOffsetCirculator<Neighborhood> nc(x == w-1
00121                                                         ? ncstartBorder
00122                                                         : ncstart);
00123             NeighborOffsetCirculator<Neighborhood> nce(x == 0
00124                                                          ? ncendBorder
00125                                                          : ncend);
00126             LabelType currentLabel = labels.nextFreeLabel();
00127             for(; nc != nce; ++nc)
00128             {
00129                 if((sa(xs) & nc.directionBit()) || (sa(xs, *nc) & nc.oppositeDirectionBit()))
00130                 {
00131                     currentLabel = labels.makeUnion(da(xd,*nc), currentLabel);
00132                 }
00133             }
00134             da.set(labels.finalizeLabel(currentLabel), xd);
00135         }
00136     }
00137 
00138     unsigned int count = labels.makeContiguous();
00139 
00140     // pass 2: assign one label to each region (tree)
00141     // so that labels form a consecutive sequence 1, 2, ...
00142     yd = upperleftd;
00143     for(y=0; y != h; ++y, ++yd.y)
00144     {
00145         DestIterator xd(yd);
00146         for(x = 0; x != w; ++x, ++xd.x)
00147         {
00148             da.set(labels[da(xd)], xd);
00149         }
00150     }
00151     return count;
00152 }
00153 
00154 template <class SrcIterator, class SrcAccessor,
00155           class DestIterator, class DestAccessor,
00156           class Neighborhood>
00157 unsigned int watershedLabeling(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00158                                pair<DestIterator, DestAccessor> dest,
00159                                Neighborhood neighborhood)
00160 {
00161     return watershedLabeling(src.first, src.second, src.third,
00162                              dest.first, dest.second, neighborhood);
00163 }
00164 
00165 
00166 template <class SrcIterator, class SrcAccessor,
00167           class DestIterator, class DestAccessor>
00168 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00169                       DestIterator upperleftd, DestAccessor da,
00170                       FourNeighborCode)
00171 {
00172     int w = lowerrights.x - upperlefts.x;
00173     int h = lowerrights.y - upperlefts.y;
00174     int x,y;
00175 
00176     SrcIterator ys(upperlefts);
00177     SrcIterator xs(ys);
00178 
00179     DestIterator yd = upperleftd;
00180 
00181     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00182     {
00183         xs = ys;
00184         DestIterator xd = yd;
00185 
00186         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00187         {
00188             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00189             typename SrcAccessor::value_type v = sa(xs);
00190             // the following choice causes minima to point
00191             // to their lowest neighbor -- would this be better???
00192             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00193             int o = 0; // means center is minimum
00194             if(atBorder == NotAtBorder)
00195             {
00196                 NeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs), cend(c);
00197                 do {
00198                     if(sa(c) <= v)
00199                     {
00200                         v = sa(c);
00201                         o = c.directionBit();
00202                     }
00203                 }
00204                 while(++c != cend);
00205             }
00206             else
00207             {
00208                 RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs, atBorder), cend(c);
00209                 do {
00210                     if(sa(c) <= v)
00211                     {
00212                         v = sa(c);
00213                         o = c.directionBit();
00214                     }
00215                 }
00216                 while(++c != cend);
00217             }
00218             da.set(o, xd);
00219         }
00220     }
00221 }
00222 
00223 template <class SrcIterator, class SrcAccessor,
00224           class DestIterator, class DestAccessor>
00225 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00226                       DestIterator upperleftd, DestAccessor da,
00227                       EightNeighborCode)
00228 {
00229     int w = lowerrights.x - upperlefts.x;
00230     int h = lowerrights.y - upperlefts.y;
00231     int x,y;
00232 
00233     SrcIterator ys(upperlefts);
00234     SrcIterator xs(ys);
00235 
00236     DestIterator yd = upperleftd;
00237 
00238     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00239     {
00240         xs = ys;
00241         DestIterator xd = yd;
00242 
00243         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00244         {
00245             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00246             typename SrcAccessor::value_type v = sa(xs);
00247             // the following choice causes minima to point
00248             // to their lowest neighbor -- would this be better???
00249             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00250             int o = 0; // means center is minimum
00251             if(atBorder == NotAtBorder)
00252             {
00253                 // handle diagonal and principal neighbors separately
00254                 // so that principal neighbors are preferred when there are
00255                 // candidates with equal strength
00256                 NeighborhoodCirculator<SrcIterator, EightNeighborCode>
00257                                       c(xs, EightNeighborCode::NorthEast);
00258                 for(int i = 0; i < 4; ++i, c += 2)
00259                 {
00260                     if(sa(c) <= v)
00261                     {
00262                         v = sa(c);
00263                         o = c.directionBit();
00264                     }
00265                 }
00266                 --c;
00267                 for(int i = 0; i < 4; ++i, c += 2)
00268                 {
00269                     if(sa(c) <= v)
00270                     {
00271                         v = sa(c);
00272                         o = c.directionBit();
00273                     }
00274                 }
00275             }
00276             else
00277             {
00278                 RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode>
00279                              c(xs, atBorder), cend(c);
00280                 do
00281                 {
00282                     if(!c.isDiagonal())
00283                         continue;
00284                     if(sa(c) <= v)
00285                     {
00286                         v = sa(c);
00287                         o = c.directionBit();
00288                     }
00289                 }
00290                 while(++c != cend);
00291                 do
00292                 {
00293                     if(c.isDiagonal())
00294                         continue;
00295                     if(sa(c) <= v)
00296                     {
00297                         v = sa(c);
00298                         o = c.directionBit();
00299                     }
00300                 }
00301                 while(++c != cend);
00302             }
00303             da.set(o, xd);
00304         }
00305     }
00306 }
00307 
00308 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00309     Region growing, watersheds, and voronoi tesselation
00310 */
00311 //@{
00312 
00313     /**\brief Options object for generateWatershedSeeds().
00314      *
00315         <b> Usage:</b>
00316 
00317         <b>\#include</b> <vigra/watersheds.hxx><br>
00318         Namespace: vigra
00319         
00320         \code
00321         IImage seeds(boundary_indicator.size());
00322         
00323         // detect all minima in 'boundary_indicator' that are below gray level 22
00324         generateWatershedSeeds(srcImageRange(boundary_indicator),
00325                                destImage(seeds),
00326                                SeedOptions().minima().threshold(22.0));
00327         \endcode
00328      */
00329 class SeedOptions
00330 {
00331 public:
00332     enum DetectMinima { LevelSets, Minima, ExtendedMinima, Unspecified };
00333     
00334     double thresh;
00335     DetectMinima mini;
00336     
00337         /**\brief Construct default options object.
00338          *
00339             Defaults are: detect minima without thresholding (i.e. all minima).
00340          */
00341     SeedOptions()
00342     : thresh(NumericTraits<double>::max()),
00343       mini(Minima)
00344     {}
00345     
00346         /** Generate seeds at minima.
00347         
00348             Default: true
00349          */
00350     SeedOptions & minima()
00351     {
00352         mini = Minima;
00353         return *this;
00354     }
00355     
00356         /** Generate seeds at minima and minimal plateaus.
00357         
00358             Default: false
00359          */
00360     SeedOptions & extendedMinima()
00361     {
00362         mini = ExtendedMinima;
00363         return *this;
00364     }
00365     
00366         /** Generate seeds as level sets.
00367         
00368             Note that you must also set a threshold to define which level set is to be used.<br>
00369             Default: false
00370          */
00371     SeedOptions & levelSets()
00372     {
00373         mini = LevelSets;
00374         return *this;
00375     }
00376     
00377         /** Generate seeds as level sets at given threshold.
00378         
00379             Equivalent to <tt>SeedOptions().levelSet().threshold(threshold)</tt><br>
00380             Default: false
00381          */
00382     SeedOptions & levelSets(double threshold)
00383     {
00384         mini = LevelSets;
00385         thresh = threshold;
00386         return *this;
00387     }
00388     
00389         /** Set threshold.
00390         
00391             The threshold will be used by both the minima and level set variants
00392             of seed generation.<br>
00393             Default: no thresholding
00394          */
00395     SeedOptions & threshold(double threshold)
00396     {
00397         thresh = threshold;
00398         return *this;
00399     }
00400     
00401         // check whether the threshold has been set for the target type T
00402     template <class T>
00403     bool thresholdIsValid() const
00404     {
00405         return thresh < double(NumericTraits<T>::max());
00406     }
00407     
00408         // indicate that this option object is invalid (for internal use in watersheds)
00409     SeedOptions & unspecified()
00410     {
00411         mini = Unspecified;
00412         return *this;
00413     }
00414 };
00415 
00416 /** \brief Generate seeds for watershed computation and seeded region growing.
00417 
00418     The source image is a boundary indicator such as the gradient magnitude
00419     or the trace of the \ref boundaryTensor(). Seeds are generally generated
00420     at locations where the boundaryness (i.e. the likelihood of the point being on the
00421     boundary) is very small. In particular, seeds can be placed by either
00422     looking for local minima (possibly including minimal plateaus) of the boundaryness,
00423     of by looking at level sets (i.e. regions where the boundaryness is below a threshold).
00424     Both methods can also be combined, so that only minima below a threshold are returned.
00425     The particular seeding strategy is specified by the <tt>options</tt> object 
00426     (see \ref SeedOptions).
00427     
00428     The pixel type of the input image must be <tt>LessThanComparable</tt>.
00429     The pixel type of the output image must be large enough to hold the labels for all seeds.
00430     (typically, you will use <tt>UInt32</tt>). The function will label seeds by consecutive integers
00431     (starting from 1) and returns the largest label it used.
00432     
00433     Pass \ref vigra::EightNeighborCode or \ref vigra::FourNeighborCode to determine the 
00434     neighborhood where pixel values are compared. 
00435     
00436     The function uses accessors.
00437 
00438     <b> Declarations:</b>
00439 
00440     pass arguments explicitly:
00441     \code
00442     namespace vigra {
00443         template <class SrcIterator, class SrcAccessor,
00444                   class DestIterator, class DestAccessor,
00445                   class Neighborhood = EightNeighborCode>
00446         unsigned int
00447         generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00448                                DestIterator upperleftd, DestAccessor da, 
00449                                Neighborhood neighborhood = EightNeighborCode(),
00450                                SeedOptions const & options = SeedOptions());
00451     }
00452     \endcode
00453 
00454     use argument objects in conjunction with \ref ArgumentObjectFactories :
00455     \code
00456     namespace vigra {
00457         template <class SrcIterator, class SrcAccessor,
00458                   class DestIterator, class DestAccessor,
00459                   class Neighborhood = EightNeighborCode>
00460         unsigned int
00461         generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00462                                pair<DestIterator, DestAccessor> dest, 
00463                                Neighborhood neighborhood = EightNeighborCode(),
00464                                SeedOptions const & options = SeedOptions());
00465     }
00466     \endcode
00467 
00468     <b> Usage:</b>
00469 
00470     <b>\#include</b> <vigra/watersheds.hxx><br>
00471     Namespace: vigra
00472 
00473     For detailed examples see watershedsRegionGrowing().
00474 */
00475 doxygen_overloaded_function(template <...> unsigned int generateWatershedSeeds)
00476 
00477 template <class SrcIterator, class SrcAccessor,
00478           class DestIterator, class DestAccessor,
00479           class Neighborhood>
00480 unsigned int
00481 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00482                        DestIterator upperleftd, DestAccessor da, 
00483                        Neighborhood neighborhood,
00484                        SeedOptions const & options = SeedOptions())
00485 {
00486     using namespace functor;
00487     typedef typename SrcAccessor::value_type SrcType;
00488     
00489     vigra_precondition(options.mini != SeedOptions::LevelSets || 
00490                        options.thresholdIsValid<SrcType>(),
00491         "generateWatershedSeeds(): SeedOptions.levelSets() must be specified with threshold.");
00492     
00493     Diff2D shape = lowerrights - upperlefts;
00494     BImage seeds(shape);
00495     
00496     if(options.mini == SeedOptions::LevelSets)
00497     {
00498         transformImage(srcIterRange(upperlefts, lowerrights, sa),
00499                        destImage(seeds),
00500                        ifThenElse(Arg1() <= Param(options.thresh), Param(1), Param(0)));
00501     }
00502     else
00503     {
00504         LocalMinmaxOptions lm_options;
00505         lm_options.neighborhood(Neighborhood::DirectionCount)
00506                   .markWith(1.0)
00507                   .allowAtBorder()
00508                   .allowPlateaus(options.mini == SeedOptions::ExtendedMinima);
00509         if(options.thresholdIsValid<SrcType>())
00510             lm_options.threshold(options.thresh);
00511             
00512         localMinima(srcIterRange(upperlefts, lowerrights, sa), destImage(seeds),
00513                     lm_options);
00514     }
00515     
00516     return labelImageWithBackground(srcImageRange(seeds), destIter(upperleftd, da), 
00517                                     Neighborhood::DirectionCount == 8, 0);
00518 }
00519 
00520 template <class SrcIterator, class SrcAccessor,
00521           class DestIterator, class DestAccessor>
00522 inline unsigned int
00523 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00524                        DestIterator upperleftd, DestAccessor da, 
00525                        SeedOptions const & options = SeedOptions())
00526 {
00527     return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da, 
00528                                    EightNeighborCode(), options);
00529 }
00530 
00531 template <class SrcIterator, class SrcAccessor,
00532           class DestIterator, class DestAccessor,
00533           class Neighborhood>
00534 inline unsigned int
00535 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00536                        pair<DestIterator, DestAccessor> dest, 
00537                        Neighborhood neighborhood,
00538                        SeedOptions const & options = SeedOptions())
00539 {
00540     return generateWatershedSeeds(src.first, src.second, src.third,
00541                                    dest.first, dest.second,    
00542                                    neighborhood, options);
00543 }
00544 
00545 template <class SrcIterator, class SrcAccessor,
00546           class DestIterator, class DestAccessor>
00547 inline unsigned int
00548 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00549                        pair<DestIterator, DestAccessor> dest, 
00550                        SeedOptions const & options = SeedOptions())
00551 {
00552     return generateWatershedSeeds(src.first, src.second, src.third,
00553                                    dest.first, dest.second,    
00554                                    EightNeighborCode(), options);
00555 }
00556 
00557 /********************************************************/
00558 /*                                                      */
00559 /*                      watersheds                      */
00560 /*                                                      */
00561 /********************************************************/
00562 
00563 /** \brief Region segmentation by means of the union-find watershed algorithm.
00564 
00565     This function implements the union-find version of the watershed algorithms
00566     as described in
00567 
00568     J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
00569     and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000
00570 
00571     The source image is a boundary indicator such as the gaussianGradientMagnitude()
00572     or the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
00573     are used as region seeds, and all other pixels are recursively assigned to the same
00574     region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
00575     \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
00576     are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
00577     The function uses accessors.
00578 
00579     Note that VIGRA provides an alternative implementation of the watershed transform via
00580     \ref watershedsRegionGrowing(). It is slower, but offers many more configuration options.
00581 
00582     <b> Declarations:</b>
00583 
00584     pass arguments explicitly:
00585     \code
00586     namespace vigra {
00587         template <class SrcIterator, class SrcAccessor,
00588                   class DestIterator, class DestAccessor,
00589                   class Neighborhood = EightNeighborCode>
00590         unsigned int
00591         watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00592                             DestIterator upperleftd, DestAccessor da,
00593                             Neighborhood neighborhood = EightNeighborCode())
00594     }
00595     \endcode
00596 
00597     use argument objects in conjunction with \ref ArgumentObjectFactories :
00598     \code
00599     namespace vigra {
00600         template <class SrcIterator, class SrcAccessor,
00601                   class DestIterator, class DestAccessor,
00602                   class Neighborhood = EightNeighborCode>
00603         unsigned int
00604         watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00605                             pair<DestIterator, DestAccessor> dest,
00606                             Neighborhood neighborhood = EightNeighborCode())
00607     }
00608     \endcode
00609 
00610     <b> Usage:</b>
00611 
00612     <b>\#include</b> <vigra/watersheds.hxx><br>
00613     Namespace: vigra
00614 
00615     Example: watersheds of the gradient magnitude.
00616 
00617     \code
00618     vigra::BImage in(w,h);
00619     ... // read input data
00620 
00621     // compute gradient magnitude as boundary indicator
00622     vigra::FImage gradMag(w, h);
00623     gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 3.0);
00624 
00625     // the pixel type of the destination image must be large enough to hold
00626     // numbers up to 'max_region_label' to prevent overflow
00627     vigra::IImage labeling(w,h);
00628     int max_region_label = watershedsUnionFind(srcImageRange(gradMag), destImage(labeling));
00629 
00630     \endcode
00631 
00632     <b> Required Interface:</b>
00633 
00634     \code
00635     SrcIterator src_upperleft, src_lowerright;
00636     DestIterator dest_upperleft;
00637 
00638     SrcAccessor src_accessor;
00639     DestAccessor dest_accessor;
00640 
00641     // compare src values
00642     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
00643 
00644     // set result
00645     int label;
00646     dest_accessor.set(label, dest_upperleft);
00647     \endcode
00648 */
00649 doxygen_overloaded_function(template <...> unsigned int watershedsUnionFind)
00650 
00651 template <class SrcIterator, class SrcAccessor,
00652           class DestIterator, class DestAccessor,
00653           class Neighborhood>
00654 unsigned int
00655 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00656                     DestIterator upperleftd, DestAccessor da, 
00657                     Neighborhood neighborhood)
00658 {
00659     SImage orientationImage(lowerrights - upperlefts);
00660 
00661     prepareWatersheds(upperlefts, lowerrights, sa,
00662                      orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
00663     return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
00664                              upperleftd, da, neighborhood);
00665 }
00666 
00667 template <class SrcIterator, class SrcAccessor,
00668           class DestIterator, class DestAccessor>
00669 inline unsigned int
00670 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00671            DestIterator upperleftd, DestAccessor da)
00672 {
00673     return watershedsUnionFind(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
00674 }
00675 
00676 template <class SrcIterator, class SrcAccessor,
00677           class DestIterator, class DestAccessor,
00678           class Neighborhood>
00679 inline unsigned int
00680 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00681            pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
00682 {
00683     return watershedsUnionFind(src.first, src.second, src.third, 
00684                                 dest.first, dest.second, neighborhood);
00685 }
00686 
00687 template <class SrcIterator, class SrcAccessor,
00688           class DestIterator, class DestAccessor>
00689 inline unsigned int
00690 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00691            pair<DestIterator, DestAccessor> dest)
00692 {
00693     return watershedsUnionFind(src.first, src.second, src.third, 
00694                                 dest.first, dest.second);
00695 }
00696 
00697 /** \brief Options object for watershedsRegionGrowing().
00698 
00699     <b> Usage:</b>
00700 
00701     see watershedsRegionGrowing() for detailed examples.
00702 */
00703 class WatershedOptions
00704 {
00705   public:
00706     double max_cost, bias;
00707     SRGType terminate;
00708     unsigned int biased_label, bucket_count;
00709     SeedOptions seed_options;
00710     
00711         /** \brief Create options object with default settings.
00712 
00713             Defaults are: perform complete grow (all pixels are assigned to regions),
00714             use standard algorithm, assume that the destination image already contains 
00715             region seeds.
00716         */
00717     WatershedOptions()
00718     : max_cost(0.0),
00719       bias(1.0),
00720       terminate(CompleteGrow),
00721       biased_label(0),
00722       bucket_count(0),
00723       seed_options(SeedOptions().unspecified())
00724     {}    
00725     
00726         /** \brief Perform complete grow.
00727 
00728             That is, all pixels are assigned to regions, without explicit contours
00729             in between.
00730             
00731             Default: true
00732         */
00733     WatershedOptions & completeGrow()
00734     {
00735         terminate = SRGType(CompleteGrow | (terminate & StopAtThreshold));
00736         return *this;
00737     }
00738     
00739         /** \brief Keep one-pixel wide contour between regions.
00740         
00741             Note that this option is unsupported by the turbo algorithm.
00742 
00743             Default: false
00744         */
00745     WatershedOptions & keepContours()
00746     {
00747         terminate = SRGType(KeepContours | (terminate & StopAtThreshold));
00748         return *this;
00749     }
00750     
00751         /** \brief Set \ref SRGType explicitly.
00752         
00753             Default: CompleteGrow
00754         */
00755     WatershedOptions & srgType(SRGType type)
00756     {
00757         terminate = type;
00758         return *this;
00759     }
00760     
00761         /** \brief Stop region growing when the boundaryness exceeds the threshold.
00762         
00763             This option may be combined with completeGrow() and keepContours().
00764         
00765             Default: no early stopping
00766         */
00767     WatershedOptions & stopAtThreshold(double threshold)
00768     {
00769         terminate = SRGType(terminate | StopAtThreshold);
00770         max_cost = threshold;
00771         return *this;
00772     }
00773     
00774         /** \brief Use a simpler, but faster region growing algorithm.
00775         
00776             The algorithm internally uses a \ref BucketQueue to determine
00777             the processing order of the pixels. This is only useful,
00778             when the input boundary indicator image contains integers
00779             in the range <tt>[0, ..., bucket_count-1]</tt>. Since
00780             these boundary indicators are typically represented as
00781             UInt8 images, the default <tt>bucket_count</tt> is 256.
00782         
00783             Default: don't use the turbo algorithm
00784         */
00785     WatershedOptions & turboAlgorithm(unsigned int bucket_count = 256)
00786     {
00787         this->bucket_count = bucket_count;
00788         return *this;
00789     }
00790     
00791         /** \brief Specify seed options.
00792         
00793             In this case, watershedsRegionGrowing() assumes that the destination
00794             image does not yet contain seeds. It will therefore call 
00795             generateWatershedSeeds() and pass on the seed options.
00796         
00797             Default: don't compute seeds (i.e. assume that destination image already
00798             contains seeds).
00799         */
00800     WatershedOptions & seedOptions(SeedOptions const & s)
00801     {
00802         seed_options = s;
00803         return *this;
00804     }
00805     
00806         /** \brief Bias the cost of the specified region by the given factor.
00807         
00808             In certain applications, one region (typically the background) should
00809             be preferred in region growing. This is most easily achieved
00810             by adjusting the assignment cost for that region as <tt>factor*cost</tt>,
00811             with a factor slightly below 1.
00812         
00813             Default: don't bias any region.
00814         */
00815     WatershedOptions & biasLabel(unsigned int label, double factor)
00816     {
00817         biased_label = label;
00818         bias = factor;
00819         return *this;
00820     }
00821 };
00822 
00823 namespace detail {
00824 
00825 template <class CostType, class LabelType>
00826 class WatershedStatistics
00827 {
00828   public:
00829   
00830     typedef SeedRgDirectValueFunctor<CostType> value_type;
00831     typedef value_type & reference;
00832     typedef value_type const & const_reference;
00833     
00834     typedef CostType  first_argument_type;
00835     typedef LabelType second_argument_type;
00836     typedef LabelType argument_type;
00837     
00838     WatershedStatistics()
00839     {}
00840 
00841     void resize(unsigned int)
00842     {}
00843 
00844     void reset()
00845     {}
00846 
00847         /** update regions statistics (do nothing in the watershed algorithm)
00848         */
00849     template <class T1, class T2>
00850     void operator()(first_argument_type const &, second_argument_type const &) 
00851     {}
00852 
00853         /** ask for maximal index (label) allowed
00854         */
00855     LabelType maxRegionLabel() const
00856         { return size() - 1; }
00857 
00858         /** ask for array size (i.e. maxRegionLabel() + 1)
00859         */
00860     LabelType size() const
00861         { return NumericTraits<LabelType>::max(); }
00862 
00863         /** read the statistics functor for a region via its label
00864         */
00865     const_reference operator[](argument_type label) const
00866         { return stats; }
00867 
00868         /** access the statistics functor for a region via its label
00869         */
00870     reference operator[](argument_type label)
00871         { return stats; }
00872 
00873     value_type stats;
00874 };
00875 
00876 template <class Value>
00877 class SeedRgBiasedValueFunctor
00878 {
00879   public:
00880     double bias;
00881 
00882         /* the functor's argument type
00883         */
00884     typedef Value argument_type;
00885 
00886         /* the functor's result type (unused, only necessary for
00887             use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
00888         */
00889     typedef Value result_type;
00890 
00891         /* the return type of the cost() function
00892         */
00893     typedef Value cost_type;
00894     
00895     SeedRgBiasedValueFunctor(double b = 1.0)
00896     : bias(b)
00897     {}
00898 
00899         /* Do nothing (since we need not update region statistics).
00900         */
00901     void operator()(argument_type const &) const {}
00902 
00903         /* Return scaled argument
00904         */
00905     cost_type cost(argument_type const & v) const
00906     {
00907         return cost_type(bias*v);
00908     }
00909 };
00910 
00911 template <class CostType, class LabelType>
00912 class BiasedWatershedStatistics
00913 {
00914   public:
00915   
00916     typedef SeedRgBiasedValueFunctor<CostType> value_type;
00917     typedef value_type & reference;
00918     typedef value_type const & const_reference;
00919     
00920     typedef CostType  first_argument_type;
00921     typedef LabelType second_argument_type;
00922     typedef LabelType argument_type;
00923     
00924     BiasedWatershedStatistics(LabelType biasedLabel, double bias)
00925     : biased_label(biasedLabel),
00926       biased_stats(bias)
00927     {}
00928 
00929     void resize(unsigned int)
00930     {}
00931 
00932     void reset()
00933     {}
00934 
00935         /** update regions statistics (do nothing in the watershed algorithm)
00936         */
00937     template <class T1, class T2>
00938     void operator()(first_argument_type const &, second_argument_type const &) 
00939     {}
00940 
00941         /** ask for maximal index (label) allowed
00942         */
00943     LabelType maxRegionLabel() const
00944         { return size() - 1; }
00945 
00946         /** ask for array size (i.e. maxRegionLabel() + 1)
00947         */
00948     LabelType size() const
00949         { return NumericTraits<LabelType>::max(); }
00950 
00951         /** read the statistics functor for a region via its label
00952         */
00953     const_reference operator[](argument_type label) const
00954     { 
00955         return (label == biased_label)
00956                     ? biased_stats
00957                     : stats; 
00958     }
00959 
00960         /** access the statistics functor for a region via its label
00961         */
00962     reference operator[](argument_type label)
00963     { 
00964         return (label == biased_label)
00965                     ? biased_stats
00966                     : stats; 
00967     }
00968 
00969     LabelType biased_label;
00970     value_type stats, biased_stats;
00971 };
00972 
00973 } // namespace detail
00974 
00975 /** \brief Region segmentation by means of a flooding-based watershed algorithm.
00976 
00977     This function implements variants of the watershed algorithm
00978     described in
00979 
00980     L. Vincent and P. Soille: "<em>Watersheds in digital spaces: An efficient algorithm
00981     based on immersion simulations</em>", IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991
00982 
00983     The source image is a boundary indicator such as the gaussianGradientMagnitude()
00984     or the trace of the \ref boundaryTensor(), and the destination is a label image
00985     designating membership of each pixel in one of the regions. Plateaus in the boundary
00986     indicator (i.e. regions of constant gray value) are handled via a Euclidean distance
00987     transform by default.
00988     
00989     By default, the destination image is assumed to hold seeds for a seeded watershed 
00990     transform. Seeds may, for example, be created by means of generateWatershedSeeds(). 
00991     Note that the seeds will be overridden with the final watershed segmentation.
00992     
00993     Alternatively, you may provide \ref SeedOptions in order to instruct 
00994     watershedsRegionGrowing() to generate its own seeds (it will call generateWatershedSeeds()
00995     internally). In that case, the destination image should be zero-initialized.
00996     
00997     You can specify the neighborhood system to be used by passing \ref FourNeighborCode 
00998     or \ref EightNeighborCode (default).
00999     
01000     Further options to be specified via \ref WatershedOptions are:
01001     
01002     <ul>
01003     <li> Whether to keep a 1-pixel-wide contour (with label 0) between regions or 
01004          perform complete grow (i.e. all pixels are assigned to a region).
01005     <li> Whether to stop growing when the boundaryness exceeds a threshold (remaining
01006          pixels keep label 0).
01007     <li> Whether to use a faster, but less powerful algorithm ("turbo algorithm"). It
01008          is faster because it orders pixels by means of a \ref BucketQueue (therefore,
01009          the boundary indicator must contain integers in the range 
01010          <tt>[0, ..., bucket_count-1]</tt>, where <tt>bucket_count</tt> is specified in
01011          the options object), it only supports complete growing (no contour between regions
01012          is possible), and it handles plateaus in a simplistic way. It also saves some
01013          memory because it allocates less temporary storage.
01014     <li> Whether one region (label) is to be preferred or discouraged by biasing its cost 
01015          with a given factor (smaller than 1 for preference, larger than 1 for discouragement).
01016     </ul>
01017 
01018     Note that VIGRA provides an alternative implementation of the watershed transform via
01019     \ref watershedsUnionFind(). 
01020 
01021     <b> Declarations:</b>
01022 
01023     pass arguments explicitly:
01024     \code
01025     namespace vigra {
01026         template <class SrcIterator, class SrcAccessor,
01027                   class DestIterator, class DestAccessor,
01028                   class Neighborhood = EightNeighborCode>
01029         unsigned int
01030         watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01031                                 DestIterator upperleftd, DestAccessor da, 
01032                                 Neighborhood neighborhood = EightNeighborCode(),
01033                                 WatershedOptions const & options = WatershedOptions());
01034 
01035         template <class SrcIterator, class SrcAccessor,
01036                   class DestIterator, class DestAccessor>
01037         unsigned int
01038         watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01039                                 DestIterator upperleftd, DestAccessor da, 
01040                                 WatershedOptions const & options = WatershedOptions());
01041     }
01042     \endcode
01043 
01044     use argument objects in conjunction with \ref ArgumentObjectFactories :
01045     \code
01046     namespace vigra {
01047         template <class SrcIterator, class SrcAccessor,
01048                   class DestIterator, class DestAccessor,
01049                   class Neighborhood = EightNeighborCode>
01050         unsigned int
01051         watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01052                                 pair<DestIterator, DestAccessor> dest, 
01053                                 Neighborhood neighborhood = EightNeighborCode(),
01054                                 WatershedOptions const & options = WatershedOptions());
01055                                 
01056         template <class SrcIterator, class SrcAccessor,
01057                   class DestIterator, class DestAccessor>
01058         unsigned int
01059         watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01060                                 pair<DestIterator, DestAccessor> dest, 
01061                                 WatershedOptions const & options = WatershedOptions());
01062     }
01063     \endcode
01064 
01065     <b> Usage:</b>
01066 
01067     <b>\#include</b> <vigra/watersheds.hxx><br>
01068     Namespace: vigra
01069 
01070     Example: watersheds of the gradient magnitude.
01071 
01072     \code
01073     vigra::BImage src(w, h);
01074     ... // read input data
01075     
01076     // compute gradient magnitude at scale 1.0 as a boundary indicator
01077     vigra::FImage gradMag(w, h);
01078     gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0);
01079 
01080     // example 1
01081     {
01082         // the pixel type of the destination image must be large enough to hold
01083         // numbers up to 'max_region_label' to prevent overflow
01084         vigra::IImage labeling(w, h);
01085         
01086         // call watershed algorithm for 4-neighborhood, leave a 1-pixel boundary between regions,
01087         // and autogenerate seeds from all gradient minima where the magnitude is below 2.0
01088         unsigned int max_region_label = 
01089               watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
01090                                       FourNeighborCode(),
01091                                       WatershedOptions().keepContours()
01092                                            .seedOptions(SeedOptions().minima().threshold(2.0)));
01093     }
01094     
01095     // example 2
01096     {
01097         vigra::IImage labeling(w, h);
01098         
01099         // compute seeds beforehand (use connected components of all pixels 
01100         // where the gradient  is below 4.0)
01101         unsigned int max_region_label = 
01102               generateWatershedSeeds(srcImageRange(gradMag), destImage(labeling),
01103                                      SeedOptions().levelSets(4.0));
01104         
01105         // quantize the gradient image to 256 gray levels
01106         vigra::BImage gradMag256(w, h);
01107         vigra::FindMinMax<float> minmax; 
01108         inspectImage(srcImageRange(gradMag), minmax); // find original range
01109         transformImage(srcImageRange(gradMag), destImage(gradMag256),
01110                        linearRangeMapping(minmax, 0, 255));
01111         
01112         // call the turbo algorithm with 256 bins, using 8-neighborhood
01113         watershedsRegionGrowing(srcImageRange(gradMag256), destImage(labeling),
01114                                 WatershedOptions().turboAlgorithm(256));
01115     }
01116     
01117     // example 3
01118     {
01119         vigra::IImage labeling(w, h);
01120         
01121         .. // get seeds from somewhere, e.g. an interactive labeling program,
01122            // make sure that label 1 corresponds to the background
01123         
01124         // bias the watershed algorithm so that the background is preferred
01125         // by reducing the cost for label 1 to 90%
01126         watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
01127                                 WatershedOptions().biasLabel(1, 0.9));
01128     }
01129     \endcode
01130 
01131     <b> Required Interface:</b>
01132 
01133     \code
01134     SrcIterator src_upperleft, src_lowerright;
01135     DestIterator dest_upperleft;
01136 
01137     SrcAccessor src_accessor;
01138     DestAccessor dest_accessor;
01139 
01140     // compare src values
01141     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
01142 
01143     // set result
01144     int label;
01145     dest_accessor.set(label, dest_upperleft);
01146     \endcode
01147 */
01148 doxygen_overloaded_function(template <...> unsigned int watershedsRegionGrowing)
01149 
01150 template <class SrcIterator, class SrcAccessor,
01151           class DestIterator, class DestAccessor,
01152           class Neighborhood>
01153 unsigned int
01154 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01155                         DestIterator upperleftd, DestAccessor da, 
01156                         Neighborhood neighborhood,
01157                         WatershedOptions const & options = WatershedOptions())
01158 {
01159     typedef typename SrcAccessor::value_type ValueType; 
01160     typedef typename DestAccessor::value_type LabelType; 
01161     
01162     unsigned int max_region_label = 0;
01163     
01164     if(options.seed_options.mini != SeedOptions::Unspecified)
01165     {
01166         // we are supposed to compute seeds
01167         max_region_label = 
01168             generateWatershedSeeds(srcIterRange(upperlefts, lowerrights, sa), 
01169                                    destIter(upperleftd, da),
01170                                    neighborhood, options.seed_options);
01171     }
01172     
01173     if(options.biased_label != 0)
01174     {
01175         // create a statistics functor for biased region growing
01176         detail::BiasedWatershedStatistics<ValueType, LabelType> 
01177                                  regionstats(options.biased_label, options.bias);
01178 
01179         // perform region growing, starting from the seeds computed above
01180         if(options.bucket_count == 0)
01181         {
01182             max_region_label = 
01183             seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01184                                 srcIter(upperleftd, da),
01185                                 destIter(upperleftd, da), 
01186                                 regionstats, options.terminate, neighborhood, options.max_cost);
01187         }
01188         else
01189         {
01190             max_region_label = 
01191             fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01192                                     destIter(upperleftd, da), 
01193                                     regionstats, options.terminate, 
01194                                     neighborhood, options.max_cost, options.bucket_count);
01195         }
01196     }
01197     else
01198     {
01199         // create a statistics functor for region growing
01200         detail::WatershedStatistics<ValueType, LabelType> regionstats;
01201 
01202         // perform region growing, starting from the seeds computed above
01203         if(options.bucket_count == 0)
01204         {
01205             max_region_label = 
01206             seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01207                                 srcIter(upperleftd, da),
01208                                 destIter(upperleftd, da), 
01209                                 regionstats, options.terminate, neighborhood, options.max_cost);
01210         }
01211         else
01212         {
01213             max_region_label = 
01214             fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01215                                     destIter(upperleftd, da), 
01216                                     regionstats, options.terminate, 
01217                                     neighborhood, options.max_cost, options.bucket_count);
01218         }
01219     }
01220     
01221     return max_region_label;
01222 }
01223 
01224 template <class SrcIterator, class SrcAccessor,
01225           class DestIterator, class DestAccessor>
01226 inline unsigned int
01227 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01228                         DestIterator upperleftd, DestAccessor da, 
01229                         WatershedOptions const & options = WatershedOptions())
01230 {
01231     return watershedsRegionGrowing(upperlefts, lowerrights, sa, upperleftd,  da,
01232                                    EightNeighborCode(), options);
01233 }
01234 
01235 template <class SrcIterator, class SrcAccessor,
01236           class DestIterator, class DestAccessor,
01237           class Neighborhood>
01238 inline unsigned int
01239 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01240                         pair<DestIterator, DestAccessor> dest, 
01241                         Neighborhood neighborhood,
01242                         WatershedOptions const & options = WatershedOptions())
01243 {
01244     return watershedsRegionGrowing(src.first, src.second, src.third,
01245                                    dest.first, dest.second,    
01246                                    neighborhood, options);
01247 }
01248 
01249 template <class SrcIterator, class SrcAccessor,
01250           class DestIterator, class DestAccessor>
01251 inline unsigned int
01252 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01253                         pair<DestIterator, DestAccessor> dest, 
01254                         WatershedOptions const & options = WatershedOptions())
01255 {
01256     return watershedsRegionGrowing(src.first, src.second, src.third,
01257                                     dest.first, dest.second,    
01258                                     EightNeighborCode(), options);
01259 }
01260 
01261 
01262 //@}
01263 
01264 } // namespace vigra
01265 
01266 #endif // VIGRA_WATERSHEDS_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)