[ 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         localMinima(srcIterRange(upperlefts, lowerrights, sa), destImage(seeds),
00505             LocalMinmaxOptions().neighborhood(Neighborhood::DirectionCount)
00506                                 .markWith(1.0)
00507                                 .threshold(options.thresh)
00508                                 .allowAtBorder()
00509                                 .allowPlateaus(options.mini == SeedOptions::ExtendedMinima));
00510     }
00511     
00512     return labelImageWithBackground(srcImageRange(seeds), destIter(upperleftd, da), 
00513                                     Neighborhood::DirectionCount == 8, 0);
00514 }
00515 
00516 template <class SrcIterator, class SrcAccessor,
00517           class DestIterator, class DestAccessor>
00518 inline unsigned int
00519 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00520                        DestIterator upperleftd, DestAccessor da, 
00521                        SeedOptions const & options = SeedOptions())
00522 {
00523     return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da, 
00524                                    EightNeighborCode(), options);
00525 }
00526 
00527 template <class SrcIterator, class SrcAccessor,
00528           class DestIterator, class DestAccessor,
00529           class Neighborhood>
00530 inline unsigned int
00531 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00532                        pair<DestIterator, DestAccessor> dest, 
00533                        Neighborhood neighborhood,
00534                        SeedOptions const & options = SeedOptions())
00535 {
00536     return generateWatershedSeeds(src.first, src.second, src.third,
00537                                    dest.first, dest.second,    
00538                                    neighborhood, options);
00539 }
00540 
00541 template <class SrcIterator, class SrcAccessor,
00542           class DestIterator, class DestAccessor>
00543 inline unsigned int
00544 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00545                        pair<DestIterator, DestAccessor> dest, 
00546                        SeedOptions const & options = SeedOptions())
00547 {
00548     return generateWatershedSeeds(src.first, src.second, src.third,
00549                                    dest.first, dest.second,    
00550                                    EightNeighborCode(), options);
00551 }
00552 
00553 /********************************************************/
00554 /*                                                      */
00555 /*                      watersheds                      */
00556 /*                                                      */
00557 /********************************************************/
00558 
00559 /** \brief Region segmentation by means of the union-find watershed algorithm.
00560 
00561     This function implements the union-find version of the watershed algorithms
00562     as described in
00563 
00564     J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
00565     and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000
00566 
00567     The source image is a boundary indicator such as the gaussianGradientMagnitude()
00568     or the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
00569     are used as region seeds, and all other pixels are recursively assigned to the same
00570     region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
00571     \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
00572     are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
00573     The function uses accessors.
00574 
00575     Note that VIGRA provides an alternative implementation of the watershed transform via
00576     \ref watershedsRegionGrowing(). It is slower, but offers many more configuration options.
00577 
00578     <b> Declarations:</b>
00579 
00580     pass arguments explicitly:
00581     \code
00582     namespace vigra {
00583         template <class SrcIterator, class SrcAccessor,
00584                   class DestIterator, class DestAccessor,
00585                   class Neighborhood = EightNeighborCode>
00586         unsigned int
00587         watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00588                             DestIterator upperleftd, DestAccessor da,
00589                             Neighborhood neighborhood = EightNeighborCode())
00590     }
00591     \endcode
00592 
00593     use argument objects in conjunction with \ref ArgumentObjectFactories :
00594     \code
00595     namespace vigra {
00596         template <class SrcIterator, class SrcAccessor,
00597                   class DestIterator, class DestAccessor,
00598                   class Neighborhood = EightNeighborCode>
00599         unsigned int
00600         watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00601                             pair<DestIterator, DestAccessor> dest,
00602                             Neighborhood neighborhood = EightNeighborCode())
00603     }
00604     \endcode
00605 
00606     <b> Usage:</b>
00607 
00608     <b>\#include</b> <vigra/watersheds.hxx><br>
00609     Namespace: vigra
00610 
00611     Example: watersheds of the gradient magnitude.
00612 
00613     \code
00614     vigra::BImage in(w,h);
00615     ... // read input data
00616 
00617     // compute gradient magnitude as boundary indicator
00618     vigra::FImage gradMag(w, h);
00619     gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 3.0);
00620 
00621     // the pixel type of the destination image must be large enough to hold
00622     // numbers up to 'max_region_label' to prevent overflow
00623     vigra::IImage labeling(w,h);
00624     int max_region_label = watershedsUnionFind(srcImageRange(gradMag), destImage(labeling));
00625 
00626     \endcode
00627 
00628     <b> Required Interface:</b>
00629 
00630     \code
00631     SrcIterator src_upperleft, src_lowerright;
00632     DestIterator dest_upperleft;
00633 
00634     SrcAccessor src_accessor;
00635     DestAccessor dest_accessor;
00636 
00637     // compare src values
00638     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
00639 
00640     // set result
00641     int label;
00642     dest_accessor.set(label, dest_upperleft);
00643     \endcode
00644 */
00645 doxygen_overloaded_function(template <...> unsigned int watershedsUnionFind)
00646 
00647 template <class SrcIterator, class SrcAccessor,
00648           class DestIterator, class DestAccessor,
00649           class Neighborhood>
00650 unsigned int
00651 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00652                     DestIterator upperleftd, DestAccessor da, 
00653                     Neighborhood neighborhood)
00654 {
00655     SImage orientationImage(lowerrights - upperlefts);
00656     SImage::traverser yo = orientationImage.upperLeft();
00657 
00658     prepareWatersheds(upperlefts, lowerrights, sa,
00659                      orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
00660     return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
00661                              upperleftd, da, neighborhood);
00662 }
00663 
00664 template <class SrcIterator, class SrcAccessor,
00665           class DestIterator, class DestAccessor>
00666 inline unsigned int
00667 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00668            DestIterator upperleftd, DestAccessor da)
00669 {
00670     return watershedsUnionFind(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
00671 }
00672 
00673 template <class SrcIterator, class SrcAccessor,
00674           class DestIterator, class DestAccessor,
00675           class Neighborhood>
00676 inline unsigned int
00677 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00678            pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
00679 {
00680     return watershedsUnionFind(src.first, src.second, src.third, 
00681                                 dest.first, dest.second, neighborhood);
00682 }
00683 
00684 template <class SrcIterator, class SrcAccessor,
00685           class DestIterator, class DestAccessor>
00686 inline unsigned int
00687 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00688            pair<DestIterator, DestAccessor> dest)
00689 {
00690     return watershedsUnionFind(src.first, src.second, src.third, 
00691                                 dest.first, dest.second);
00692 }
00693 
00694 /** \brief Options object for watershedsRegionGrowing().
00695 
00696     <b> Usage:</b>
00697 
00698     see watershedsRegionGrowing() for detailed examples.
00699 */
00700 class WatershedOptions
00701 {
00702   public:
00703     double max_cost, bias;
00704     SRGType terminate;
00705     unsigned int biased_label, bucket_count;
00706     SeedOptions seed_options;
00707     
00708         /** \brief Create options object with default settings.
00709 
00710             Defaults are: perform complete grow (all pixels are assigned to regions),
00711             use standard algorithm, assume that the destination image already contains 
00712             region seeds.
00713         */
00714     WatershedOptions()
00715     : max_cost(0.0),
00716       bias(1.0),
00717       terminate(CompleteGrow),
00718       biased_label(0),
00719       bucket_count(0),
00720       seed_options(SeedOptions().unspecified())
00721     {}    
00722     
00723         /** \brief Perform complete grow.
00724 
00725             That is, all pixels are assigned to regions, without explicit contours
00726             in between.
00727             
00728             Default: true
00729         */
00730     WatershedOptions & completeGrow()
00731     {
00732         terminate = SRGType(CompleteGrow | (terminate & StopAtThreshold));
00733         return *this;
00734     }
00735     
00736         /** \brief Keep one-pixel wide contour between regions.
00737         
00738             Note that this option is unsupported by the turbo algorithm.
00739 
00740             Default: false
00741         */
00742     WatershedOptions & keepContours()
00743     {
00744         terminate = SRGType(KeepContours | (terminate & StopAtThreshold));
00745         return *this;
00746     }
00747     
00748         /** \brief Set \ref SRGType explicitly.
00749         
00750             Default: CompleteGrow
00751         */
00752     WatershedOptions & srgType(SRGType type)
00753     {
00754         terminate = type;
00755         return *this;
00756     }
00757     
00758         /** \brief Stop region growing when the boundaryness exceeds the threshold.
00759         
00760             This option may be combined with completeGrow() and keepContours().
00761         
00762             Default: no early stopping
00763         */
00764     WatershedOptions & stopAtThreshold(double threshold)
00765     {
00766         terminate = SRGType(terminate | StopAtThreshold);
00767         max_cost = threshold;
00768         return *this;
00769     }
00770     
00771         /** \brief Use a simpler, but faster region growing algorithm.
00772         
00773             The algorithm internally uses a \ref BucketQueue to determine
00774             the processing order of the pixels. This is only useful,
00775             when the input boundary indicator image contains integers
00776             in the range <tt>[0, ..., bucket_count-1]</tt>. Since
00777             these boundary indicators are typically represented as
00778             UInt8 images, the default <tt>bucket_count</tt> is 256.
00779         
00780             Default: don't use the turbo algorithm
00781         */
00782     WatershedOptions & turboAlgorithm(unsigned int bucket_count = 256)
00783     {
00784         this->bucket_count = bucket_count;
00785         return *this;
00786     }
00787     
00788         /** \brief Specify seed options.
00789         
00790             In this case, watershedsRegionGrowing() assumes that the destination
00791             image does not yet contain seeds. It will therefore call 
00792             generateWatershedSeeds() and pass on the seed options.
00793         
00794             Default: don't compute seeds (i.e. assume that destination image already
00795             contains seeds).
00796         */
00797     WatershedOptions & seedOptions(SeedOptions const & s)
00798     {
00799         seed_options = s;
00800         return *this;
00801     }
00802     
00803         /** \brief Bias the cost of the specified region by the given factor.
00804         
00805             In certain applications, one region (typically the background) should
00806             be preferred in region growing. This is most easily achieved
00807             by adjusting the assignment cost for that region as <tt>factor*cost</tt>,
00808             with a factor slightly below 1.
00809         
00810             Default: don't bias any region.
00811         */
00812     WatershedOptions & biasLabel(unsigned int label, double factor)
00813     {
00814         biased_label = label;
00815         bias = factor;
00816         return *this;
00817     }
00818 };
00819 
00820 namespace detail {
00821 
00822 template <class CostType, class LabelType>
00823 class WatershedStatistics
00824 {
00825   public:
00826   
00827     typedef SeedRgDirectValueFunctor<CostType> value_type;
00828     typedef value_type & reference;
00829     typedef value_type const & const_reference;
00830     
00831     typedef CostType  first_argument_type;
00832     typedef LabelType second_argument_type;
00833     typedef LabelType argument_type;
00834     
00835     WatershedStatistics()
00836     {}
00837 
00838     void resize(unsigned int)
00839     {}
00840 
00841     void reset()
00842     {}
00843 
00844         /** update regions statistics (do nothing in the watershed algorithm)
00845         */
00846     template <class T1, class T2>
00847     void operator()(first_argument_type const &, second_argument_type const &) 
00848     {}
00849 
00850         /** ask for maximal index (label) allowed
00851         */
00852     LabelType maxRegionLabel() const
00853         { return size() - 1; }
00854 
00855         /** ask for array size (i.e. maxRegionLabel() + 1)
00856         */
00857     LabelType size() const
00858         { return NumericTraits<LabelType>::max(); }
00859 
00860         /** read the statistics functor for a region via its label
00861         */
00862     const_reference operator[](argument_type label) const
00863         { return stats; }
00864 
00865         /** access the statistics functor for a region via its label
00866         */
00867     reference operator[](argument_type label)
00868         { return stats; }
00869 
00870     value_type stats;
00871 };
00872 
00873 template <class Value>
00874 class SeedRgBiasedValueFunctor
00875 {
00876   public:
00877     double bias;
00878 
00879         /* the functor's argument type
00880         */
00881     typedef Value argument_type;
00882 
00883         /* the functor's result type (unused, only necessary for
00884             use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
00885         */
00886     typedef Value result_type;
00887 
00888         /* the return type of the cost() function
00889         */
00890     typedef Value cost_type;
00891     
00892     SeedRgBiasedValueFunctor(double b = 1.0)
00893     : bias(b)
00894     {}
00895 
00896         /* Do nothing (since we need not update region statistics).
00897         */
00898     void operator()(argument_type const &) const {}
00899 
00900         /* Return scaled argument
00901         */
00902     cost_type cost(argument_type const & v) const
00903     {
00904         return cost_type(bias*v);
00905     }
00906 };
00907 
00908 template <class CostType, class LabelType>
00909 class BiasedWatershedStatistics
00910 {
00911   public:
00912   
00913     typedef SeedRgBiasedValueFunctor<CostType> value_type;
00914     typedef value_type & reference;
00915     typedef value_type const & const_reference;
00916     
00917     typedef CostType  first_argument_type;
00918     typedef LabelType second_argument_type;
00919     typedef LabelType argument_type;
00920     
00921     BiasedWatershedStatistics(LabelType biasedLabel, double bias)
00922     : biased_label(biasedLabel),
00923       biased_stats(bias)
00924     {}
00925 
00926     void resize(unsigned int)
00927     {}
00928 
00929     void reset()
00930     {}
00931 
00932         /** update regions statistics (do nothing in the watershed algorithm)
00933         */
00934     template <class T1, class T2>
00935     void operator()(first_argument_type const &, second_argument_type const &) 
00936     {}
00937 
00938         /** ask for maximal index (label) allowed
00939         */
00940     LabelType maxRegionLabel() const
00941         { return size() - 1; }
00942 
00943         /** ask for array size (i.e. maxRegionLabel() + 1)
00944         */
00945     LabelType size() const
00946         { return NumericTraits<LabelType>::max(); }
00947 
00948         /** read the statistics functor for a region via its label
00949         */
00950     const_reference operator[](argument_type label) const
00951     { 
00952         return (label == biased_label)
00953                     ? biased_stats
00954                     : stats; 
00955     }
00956 
00957         /** access the statistics functor for a region via its label
00958         */
00959     reference operator[](argument_type label)
00960     { 
00961         return (label == biased_label)
00962                     ? biased_stats
00963                     : stats; 
00964     }
00965 
00966     LabelType biased_label;
00967     value_type stats, biased_stats;
00968 };
00969 
00970 } // namespace detail
00971 
00972 /** \brief Region segmentation by means of a flooding-based watershed algorithm.
00973 
00974     This function implements variants of the watershed algorithm
00975     described in
00976 
00977     L. Vincent and P. Soille: "<em>Watersheds in digital spaces: An efficient algorithm
00978     based on immersion simulations</em>", IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991
00979 
00980     The source image is a boundary indicator such as the gaussianGradientMagnitude()
00981     or the trace of the \ref boundaryTensor(), and the destination is a label image
00982     designating membership of each pixel in one of the regions. Plateaus in the boundary
00983     indicator (i.e. regions of constant gray value) are handled via a Euclidean distance
00984     transform by default.
00985     
00986     By default, the destination image is assumed to hold seeds for a seeded watershed 
00987     transform. Seeds may, for example, be created by means of generateWatershedSeeds(). 
00988     Note that the seeds will be overridden with the final watershed segmentation.
00989     
00990     Alternatively, you may provide \ref SeedOptions in order to instruct 
00991     watershedsRegionGrowing() to generate its own seeds (it will call generateWatershedSeeds()
00992     internally). In that case, the destination image should be zero-initialized.
00993     
00994     You can specify the neighborhood system to be used by passing \ref FourNeighborCode 
00995     or \ref EightNeighborCode (default).
00996     
00997     Further options to be specified via \ref WatershedOptions are:
00998     
00999     <ul>
01000     <li> Whether to keep a 1-pixel-wide contour (with label 0) between regions or 
01001          perform complete grow (i.e. all pixels are assigned to a region).
01002     <li> Whether to stop growing when the boundaryness exceeds a threshold (remaining
01003          pixels keep label 0).
01004     <li> Whether to use a faster, but less powerful algorithm ("turbo algorithm"). It
01005          is faster because it orders pixels by means of a \ref BucketQueue (therefore,
01006          the boundary indicator must contain integers in the range 
01007          <tt>[0, ..., bucket_count-1]</tt>, where <tt>bucket_count</tt> is specified in
01008          the options object), it only supports complete growing (no contour between regions
01009          is possible), and it handles plateaus in a simplistic way. It also saves some
01010          memory because it allocates less temporary storage.
01011     <li> Whether one region (label) is to be preferred or discouraged by biasing its cost 
01012          with a given factor (smaller than 1 for preference, larger than 1 for discouragement).
01013     </ul>
01014 
01015     Note that VIGRA provides an alternative implementation of the watershed transform via
01016     \ref watershedsUnionFind(). 
01017 
01018     <b> Declarations:</b>
01019 
01020     pass arguments explicitly:
01021     \code
01022     namespace vigra {
01023         template <class SrcIterator, class SrcAccessor,
01024                   class DestIterator, class DestAccessor,
01025                   class Neighborhood = EightNeighborCode>
01026         unsigned int
01027         watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01028                                 DestIterator upperleftd, DestAccessor da, 
01029                                 Neighborhood neighborhood = EightNeighborCode(),
01030                                 WatershedOptions const & options = WatershedOptions());
01031 
01032         template <class SrcIterator, class SrcAccessor,
01033                   class DestIterator, class DestAccessor>
01034         unsigned int
01035         watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01036                                 DestIterator upperleftd, DestAccessor da, 
01037                                 WatershedOptions const & options = WatershedOptions());
01038     }
01039     \endcode
01040 
01041     use argument objects in conjunction with \ref ArgumentObjectFactories :
01042     \code
01043     namespace vigra {
01044         template <class SrcIterator, class SrcAccessor,
01045                   class DestIterator, class DestAccessor,
01046                   class Neighborhood = EightNeighborCode>
01047         unsigned int
01048         watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01049                                 pair<DestIterator, DestAccessor> dest, 
01050                                 Neighborhood neighborhood = EightNeighborCode(),
01051                                 WatershedOptions const & options = WatershedOptions());
01052                                 
01053         template <class SrcIterator, class SrcAccessor,
01054                   class DestIterator, class DestAccessor>
01055         unsigned int
01056         watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01057                                 pair<DestIterator, DestAccessor> dest, 
01058                                 WatershedOptions const & options = WatershedOptions());
01059     }
01060     \endcode
01061 
01062     <b> Usage:</b>
01063 
01064     <b>\#include</b> <vigra/watersheds.hxx><br>
01065     Namespace: vigra
01066 
01067     Example: watersheds of the gradient magnitude.
01068 
01069     \code
01070     vigra::BImage src(w, h);
01071     ... // read input data
01072     
01073     // compute gradient magnitude at scale 1.0 as a boundary indicator
01074     vigra::FImage gradMag(w, h);
01075     gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0);
01076 
01077     // example 1
01078     {
01079         // the pixel type of the destination image must be large enough to hold
01080         // numbers up to 'max_region_label' to prevent overflow
01081         vigra::IImage labeling(w, h);
01082         
01083         // call watershed algorithm for 4-neighborhood, leave a 1-pixel boundary between regions,
01084         // and autogenerate seeds from all gradient minima where the magnitude is below 2.0
01085         unsigned int max_region_label = 
01086               watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
01087                                       FourNeighborCode(),
01088                                       WatershedOptions().keepContours()
01089                                            .seedOptions(SeedOptions().minima().threshold(2.0)));
01090     }
01091     
01092     // example 2
01093     {
01094         vigra::IImage labeling(w, h);
01095         
01096         // compute seeds beforehand (use connected components of all pixels 
01097         // where the gradient  is below 4.0)
01098         unsigned int max_region_label = 
01099               generateWatershedSeeds(srcImageRange(gradMag), destImage(labeling),
01100                                      SeedOptions().levelSets(4.0));
01101         
01102         // quantize the gradient image to 256 gray levels
01103         vigra::BImage gradMag256(w, h);
01104         vigra::FindMinMax<float> minmax; 
01105         inspectImage(srcImageRange(gradMag), minmax); // find original range
01106         transformImage(srcImageRange(gradMag), destImage(gradMag256),
01107                        linearRangeMapping(minmax, 0, 255));
01108         
01109         // call the turbo algorithm with 256 bins, using 8-neighborhood
01110         watershedsRegionGrowing(srcImageRange(gradMag256), destImage(labeling),
01111                                 WatershedOptions().turboAlgorithm(256));
01112     }
01113     
01114     // example 3
01115     {
01116         vigra::IImage labeling(w, h);
01117         
01118         .. // get seeds from somewhere, e.g. an interactive labeling program,
01119            // make sure that label 1 corresponds to the background
01120         
01121         // bias the watershed algorithm so that the background is preferred
01122         // by reducing the cost for label 1 to 90%
01123         watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
01124                                 WatershedOptions().biasLabel(1, 0.9));
01125     }
01126     \endcode
01127 
01128     <b> Required Interface:</b>
01129 
01130     \code
01131     SrcIterator src_upperleft, src_lowerright;
01132     DestIterator dest_upperleft;
01133 
01134     SrcAccessor src_accessor;
01135     DestAccessor dest_accessor;
01136 
01137     // compare src values
01138     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
01139 
01140     // set result
01141     int label;
01142     dest_accessor.set(label, dest_upperleft);
01143     \endcode
01144 */
01145 doxygen_overloaded_function(template <...> unsigned int watershedsRegionGrowing)
01146 
01147 template <class SrcIterator, class SrcAccessor,
01148           class DestIterator, class DestAccessor,
01149           class Neighborhood>
01150 unsigned int
01151 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01152                         DestIterator upperleftd, DestAccessor da, 
01153                         Neighborhood neighborhood,
01154                         WatershedOptions const & options = WatershedOptions())
01155 {
01156     typedef typename SrcAccessor::value_type ValueType; 
01157     typedef typename DestAccessor::value_type LabelType; 
01158     
01159     unsigned int max_region_label = 0;
01160     
01161     if(options.seed_options.mini != SeedOptions::Unspecified)
01162     {
01163         // we are supposed to compute seeds
01164         max_region_label = 
01165             generateWatershedSeeds(srcIterRange(upperlefts, lowerrights, sa), 
01166                                    destIter(upperleftd, da),
01167                                    neighborhood, options.seed_options);
01168     }
01169     
01170     if(options.biased_label != 0)
01171     {
01172         // create a statistics functor for biased region growing
01173         detail::BiasedWatershedStatistics<ValueType, LabelType> 
01174                                  regionstats(options.biased_label, options.bias);
01175 
01176         // perform region growing, starting from the seeds computed above
01177         if(options.bucket_count == 0)
01178         {
01179             max_region_label = 
01180             seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01181                                 srcIter(upperleftd, da),
01182                                 destIter(upperleftd, da), 
01183                                 regionstats, options.terminate, neighborhood, options.max_cost);
01184         }
01185         else
01186         {
01187             max_region_label = 
01188             fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01189                                     destIter(upperleftd, da), 
01190                                     regionstats, options.terminate, 
01191                                     neighborhood, options.max_cost, options.bucket_count);
01192         }
01193     }
01194     else
01195     {
01196         // create a statistics functor for region growing
01197         detail::WatershedStatistics<ValueType, LabelType> regionstats;
01198 
01199         // perform region growing, starting from the seeds computed above
01200         if(options.bucket_count == 0)
01201         {
01202             max_region_label = 
01203             seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01204                                 srcIter(upperleftd, da),
01205                                 destIter(upperleftd, da), 
01206                                 regionstats, options.terminate, neighborhood, options.max_cost);
01207         }
01208         else
01209         {
01210             max_region_label = 
01211             fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
01212                                     destIter(upperleftd, da), 
01213                                     regionstats, options.terminate, 
01214                                     neighborhood, options.max_cost, options.bucket_count);
01215         }
01216     }
01217     
01218     return max_region_label;
01219 }
01220 
01221 template <class SrcIterator, class SrcAccessor,
01222           class DestIterator, class DestAccessor>
01223 inline unsigned int
01224 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
01225                         DestIterator upperleftd, DestAccessor da, 
01226                         WatershedOptions const & options = WatershedOptions())
01227 {
01228     return watershedsRegionGrowing(upperlefts, lowerrights, sa, upperleftd,  da,
01229                                    EightNeighborCode(), options);
01230 }
01231 
01232 template <class SrcIterator, class SrcAccessor,
01233           class DestIterator, class DestAccessor,
01234           class Neighborhood>
01235 inline unsigned int
01236 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01237                         pair<DestIterator, DestAccessor> dest, 
01238                         Neighborhood neighborhood,
01239                         WatershedOptions const & options = WatershedOptions())
01240 {
01241     return watershedsRegionGrowing(src.first, src.second, src.third,
01242                                    dest.first, dest.second,    
01243                                    neighborhood, options);
01244 }
01245 
01246 template <class SrcIterator, class SrcAccessor,
01247           class DestIterator, class DestAccessor>
01248 inline unsigned int
01249 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01250                         pair<DestIterator, DestAccessor> dest, 
01251                         WatershedOptions const & options = WatershedOptions())
01252 {
01253     return watershedsRegionGrowing(src.first, src.second, src.third,
01254                                     dest.first, dest.second,    
01255                                     EightNeighborCode(), options);
01256 }
01257 
01258 
01259 //@}
01260 
01261 } // namespace vigra
01262 
01263 #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.8.0 (20 Sep 2011)