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

vigra/watersheds.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2005 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    The VIGRA Website is                                              */
00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00034 /*                                                                      */
00035 /************************************************************************/
00036 
00037 #ifndef VIGRA_WATERSHEDS_HXX
00038 #define VIGRA_WATERSHEDS_HXX
00039 
00040 #include <functional>
00041 #include "mathutil.hxx"
00042 #include "stdimage.hxx"
00043 #include "pixelneighborhood.hxx"
00044 
00045 namespace vigra {
00046 
00047 template <class SrcIterator, class SrcAccessor,
00048           class DestIterator, class DestAccessor,
00049           class Neighborhood>
00050 unsigned int watershedLabeling(SrcIterator upperlefts,
00051                         SrcIterator lowerrights, SrcAccessor sa,
00052                         DestIterator upperleftd, DestAccessor da,
00053                         Neighborhood neighborhood)
00054 {
00055     int w = lowerrights.x - upperlefts.x;
00056     int h = lowerrights.y - upperlefts.y;
00057     int x,y,i;
00058 
00059     SrcIterator ys(upperlefts);
00060     SrcIterator xs(ys);
00061 
00062     // temporary image to store region labels
00063     typedef IImage LabelImage;
00064     typedef LabelImage::traverser LabelTraverser;
00065     
00066     LabelImage labelimage(w, h);
00067 
00068     LabelTraverser yt = labelimage.upperLeft();
00069     LabelTraverser xt(yt);
00070 
00071     // Kovalevsky's clever idea to use
00072     // image iterator and scan order iterator simultaneously
00073     LabelImage::ScanOrderIterator label = labelimage.begin();
00074     
00075     // initialize the neighborhood circulators
00076     NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst);
00077     NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North);
00078     NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast);
00079     ++ncend;
00080     NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North);
00081     ++ncendBorder;
00082 
00083     // pass 1: scan image from upper left to lower right
00084     // to find connected components
00085 
00086     // Each component will be represented by a tree of pixels. Each
00087     // pixel contains the scan order address of its parent in the
00088     // tree.  In order for pass 2 to work correctly, the parent must
00089     // always have a smaller scan order address than the child.
00090     // Therefore, we can merge trees only at their roots, because the
00091     // root of the combined tree must have the smallest scan order
00092     // address among all the tree's pixels/ nodes.  The root of each
00093     // tree is distinguished by pointing to itself (it contains its
00094     // own scan order address). This condition is enforced whenever a
00095     // new region is found or two regions are merged
00096     xs = ys;
00097     xt = yt;
00098     
00099     *xt = 0;
00100     
00101     ++xs.x;
00102     ++xt.x;
00103     for(x = 1; x != w; ++x, ++xs.x, ++xt.x)
00104     {
00105         if((*xs & Neighborhood::directionBit(Neighborhood::West)) ||
00106            (xs[Neighborhood::west()] & Neighborhood::directionBit(Neighborhood::East)))
00107         {
00108             *xt = xt[Neighborhood::west()];
00109         }
00110         else
00111         {
00112             *xt = x;
00113         }
00114     }
00115     
00116     ++ys.y;
00117     ++yt.y;
00118     for(y = 1; y != h; ++y, ++ys.y, ++yt.y)
00119     {
00120         xs = ys;
00121         xt = yt;
00122 
00123         for(x = 0; x != w; ++x, ++xs.x, ++xt.x)
00124         {
00125             NeighborOffsetCirculator<Neighborhood> nc(x == w-1
00126                                                         ? ncstartBorder
00127                                                         : ncstart);
00128             NeighborOffsetCirculator<Neighborhood> nce(x == 0 
00129                                                          ? ncendBorder 
00130                                                          : ncend);
00131             *xt = x + w*y; // default: new region            
00132             for(; nc != nce; ++nc)
00133             {
00134                 if((*xs & nc.directionBit()) || (xs[*nc] & nc.oppositeDirectionBit()))
00135                 {
00136                     int neighborLabel = xt[*nc];
00137                     // find the root label of a label tree
00138                     while(neighborLabel != label[neighborLabel])
00139                     {
00140                         neighborLabel = label[neighborLabel];
00141                     }
00142                     if(neighborLabel < *xt) // always keep the smallest among the possible labels
00143                     {
00144                         label[*xt] = neighborLabel;
00145                         *xt = neighborLabel;
00146                     }
00147                     else
00148                     {
00149                         label[neighborLabel] = *xt;
00150                     }
00151                 }
00152             }
00153         }
00154     }
00155 
00156     // pass 2: assign one label to each region (tree)
00157     // so that labels form a consecutive sequence 1, 2, ...
00158     DestIterator yd(upperleftd);
00159 
00160     unsigned int count = 0;
00161     i = 0;
00162     for(y=0; y != h; ++y, ++yd.y)
00163     {
00164         DestIterator xd(yd);
00165         for(x = 0; x != w; ++x, ++xd.x, ++i)
00166         {
00167             if(label[i] == i)
00168             {
00169                 label[i] = ++count;
00170             }
00171             else
00172             {
00173                 label[i] = label[label[i]];
00174             }
00175             da.set(label[i], xd);
00176         }
00177     }
00178     return count;
00179 }
00180 
00181 template <class SrcIterator, class SrcAccessor,
00182           class DestIterator, class DestAccessor>
00183 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00184                       DestIterator upperleftd, DestAccessor da, 
00185                       FourNeighborCode)
00186 {
00187     int w = lowerrights.x - upperlefts.x;
00188     int h = lowerrights.y - upperlefts.y;
00189     int x,y;
00190 
00191     SrcIterator ys(upperlefts);
00192     SrcIterator xs(ys);
00193 
00194     DestIterator yd = upperleftd;
00195 
00196     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00197     {
00198         xs = ys;
00199         DestIterator xd = yd;
00200 
00201         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00202         {
00203             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00204             typename SrcAccessor::value_type v = sa(xs);
00205             // the following choice causes minima to point
00206             // to their lowest neighbor -- would this be better???
00207             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00208             int o = 0; // means center is minimum
00209             if(atBorder == NotAtBorder)
00210             {
00211                 NeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs), cend(c);
00212                 do {
00213                     if(sa(c) <= v)
00214                     {
00215                         v = sa(c);
00216                         o = c.directionBit();
00217                     }
00218                 }
00219                 while(++c != cend);
00220             }
00221             else
00222             {
00223                 RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs, atBorder), cend(c);
00224                 do {
00225                     if(sa(c) <= v)
00226                     {
00227                         v = sa(c);
00228                         o = c.directionBit();
00229                     }
00230                 }
00231                 while(++c != cend);
00232             }
00233             da.set(o, xd);
00234         }
00235     }
00236 }
00237 
00238 template <class SrcIterator, class SrcAccessor,
00239           class DestIterator, class DestAccessor>
00240 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00241                       DestIterator upperleftd, DestAccessor da, 
00242                       EightNeighborCode)
00243 {
00244     int w = lowerrights.x - upperlefts.x;
00245     int h = lowerrights.y - upperlefts.y;
00246     int x,y;
00247 
00248     SrcIterator ys(upperlefts);
00249     SrcIterator xs(ys);
00250 
00251     DestIterator yd = upperleftd;
00252 
00253     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00254     {
00255         xs = ys;
00256         DestIterator xd = yd;
00257 
00258         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00259         {
00260             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00261             typename SrcAccessor::value_type v = sa(xs);
00262             // the following choice causes minima to point
00263             // to their lowest neighbor -- would this be better???
00264             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00265             int o = 0; // means center is minimum
00266             if(atBorder == NotAtBorder)
00267             {
00268                 // handle diagonal and principal neighbors separately
00269                 // so that principal neighbors are preferred when there are 
00270                 // candidates with equal strength
00271                 NeighborhoodCirculator<SrcIterator, EightNeighborCode>  
00272                                       c(xs, EightNeighborCode::NorthEast);
00273                 for(int i = 0; i < 4; ++i, c += 2)
00274                 {
00275                     if(sa(c) <= v)
00276                     {
00277                         v = sa(c);
00278                         o = c.directionBit();
00279                     }
00280                 }
00281                 --c;
00282                 for(int i = 0; i < 4; ++i, c += 2)
00283                 {
00284                     if(sa(c) <= v)
00285                     {
00286                         v = sa(c);
00287                         o = c.directionBit();
00288                     }
00289                 }
00290             }
00291             else
00292             {
00293                 RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode>  
00294                              c(xs, atBorder), cend(c);
00295                 do 
00296                 {
00297                     if(!c.isDiagonal())
00298                         continue;
00299                     if(sa(c) <= v)
00300                     {
00301                         v = sa(c);
00302                         o = c.directionBit();
00303                     }
00304                 }
00305                 while(++c != cend);
00306                 do 
00307                 {
00308                     if(c.isDiagonal())
00309                         continue;
00310                     if(sa(c) <= v)
00311                     {
00312                         v = sa(c);
00313                         o = c.directionBit();
00314                     }
00315                 }
00316                 while(++c != cend);
00317             }
00318             da.set(o, xd);
00319         }
00320     }
00321 }
00322 
00323 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00324     Region growing, watersheds, and voronoi tesselation
00325 */
00326 //@{
00327 
00328 /********************************************************/
00329 /*                                                      */
00330 /*                      watersheds                      */
00331 /*                                                      */
00332 /********************************************************/
00333 
00334 /** \brief Region Segmentation by means of the watershed algorithm.
00335 
00336     This function implements the union-find version of the watershed algorithms
00337     as described in
00338 
00339     J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
00340     and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000
00341 
00342     The source image is a boundary indicator such as the gradient magnitude
00343     of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
00344     are used as region seeds, and all other pixels are recursively assigned to the same 
00345     region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or 
00346     \ref vigra::FourNeighborCode to determine the neighborhood where pixel values 
00347     are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
00348     The function uses accessors. 
00349     
00350     Note that VIGRA provides an alternative implementaion of the watershed transform via
00351     \ref seededRegionGrowing(). It is slower, but handles plateaus better 
00352     and allows to keep a one pixel wide boundary between regions.
00353     
00354     <b> Declarations:</b>
00355 
00356     pass arguments explicitly:
00357     \code
00358     namespace vigra {
00359         template <class SrcIterator, class SrcAccessor,
00360                   class DestIterator, class DestAccessor,
00361                   class Neighborhood = EightNeighborCode>
00362         unsigned int 
00363         watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00364                    DestIterator upperleftd, DestAccessor da, 
00365                    Neighborhood neighborhood = EightNeighborCode())
00366     }
00367     \endcode
00368 
00369     use argument objects in conjunction with \ref ArgumentObjectFactories :
00370     \code
00371     namespace vigra {
00372         template <class SrcIterator, class SrcAccessor,
00373                   class DestIterator, class DestAccessor,
00374                   class Neighborhood = EightNeighborCode>
00375         unsigned int 
00376         watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00377                    pair<DestIterator, DestAccessor> dest, 
00378                    Neighborhood neighborhood = EightNeighborCode())
00379     }
00380     \endcode
00381 
00382     <b> Usage:</b>
00383 
00384     <b>\#include</b> <<a href="watersheds_8hxx-source.html">vigra/watersheds.hxx</a>><br>
00385     Namespace: vigra
00386 
00387     Example: watersheds of the gradient magnitude.
00388 
00389     \code
00390     vigra::BImage in(w,h);
00391     ... // read input data
00392     
00393     vigra::FImage gradx(x,y), grady(x,y), gradMag(x,y);
00394     gaussianGradient(srcImageRange(src), destImage(gradx), destImage(grady), 3.0);
00395     combineTwoImages(srcImageRange(gradx), srcImage(grady), destImage(gradMag),
00396                      vigra::MagnitudeFunctor<float>());
00397     
00398     // the pixel type of the destination image must be large enough to hold
00399     // numbers up to 'max_region_label' to prevent overflow
00400     vigra::IImage labeling(x,y);
00401     int max_region_label = watersheds(srcImageRange(gradMag), destImage(labeling));
00402 
00403     \endcode
00404 
00405     <b> Required Interface:</b>
00406 
00407     \code
00408     SrcImageIterator src_upperleft, src_lowerright;
00409     DestImageIterator dest_upperleft;
00410 
00411     SrcAccessor src_accessor;
00412     DestAccessor dest_accessor;
00413     
00414     // compare src values
00415     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
00416 
00417     // set result
00418     int label;
00419     dest_accessor.set(label, dest_upperleft);
00420     \endcode
00421 */
00422 doxygen_overloaded_function(template <...> unsigned int watersheds)
00423 
00424 template <class SrcIterator, class SrcAccessor,
00425           class DestIterator, class DestAccessor,
00426           class Neighborhood>
00427 unsigned int 
00428 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00429            DestIterator upperleftd, DestAccessor da, Neighborhood neighborhood)
00430 {
00431     SImage orientationImage(lowerrights - upperlefts);
00432     SImage::traverser yo = orientationImage.upperLeft();
00433 
00434     prepareWatersheds(upperlefts, lowerrights, sa, 
00435                      orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
00436     return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
00437                              upperleftd, da, neighborhood);
00438 }
00439 
00440 template <class SrcIterator, class SrcAccessor,
00441           class DestIterator, class DestAccessor>
00442 inline unsigned int 
00443 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00444            DestIterator upperleftd, DestAccessor da)
00445 {
00446     return watersheds(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
00447 }
00448 
00449 template <class SrcIterator, class SrcAccessor,
00450           class DestIterator, class DestAccessor,
00451           class Neighborhood>
00452 inline unsigned int 
00453 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00454            pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
00455 {
00456     return watersheds(src.first, src.second, src.third, dest.first, dest.second, neighborhood);
00457 }
00458 
00459 template <class SrcIterator, class SrcAccessor,
00460           class DestIterator, class DestAccessor>
00461 inline unsigned int 
00462 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00463            pair<DestIterator, DestAccessor> dest)
00464 {
00465     return watersheds(src.first, src.second, src.third, dest.first, dest.second);
00466 }
00467 
00468 //@}
00469 
00470 } // namespace vigra
00471 
00472 #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.6.0 (5 Nov 2009)