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

vigra/watersheds.hxx
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 "union_find.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)
00054 {
00055     typedef typename DestAccessor::value_type LabelType;
00056 
00057     int w = lowerrights.x - upperlefts.x;
00058     int h = lowerrights.y - upperlefts.y;
00059     int x,y;
00060 
00061     SrcIterator ys(upperlefts);
00062     SrcIterator xs(ys);
00063     DestIterator yd(upperleftd);
00064     DestIterator xd(yd);
00065 
00066     // temporary image to store region labels
00067     detail::UnionFindArray<LabelType> labels;
00068 
00069     // initialize the neighborhood circulators
00070     NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst);
00071     NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North);
00072     NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast);
00073     ++ncend;
00074     NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North);
00075     ++ncendBorder;
00076 
00077     // pass 1: scan image from upper left to lower right
00078     // to find connected components
00079 
00080     // Each component will be represented by a tree of pixels. Each
00081     // pixel contains the scan order address of its parent in the
00082     // tree.  In order for pass 2 to work correctly, the parent must
00083     // always have a smaller scan order address than the child.
00084     // Therefore, we can merge trees only at their roots, because the
00085     // root of the combined tree must have the smallest scan order
00086     // address among all the tree's pixels/ nodes.  The root of each
00087     // tree is distinguished by pointing to itself (it contains its
00088     // own scan order address). This condition is enforced whenever a
00089     // new region is found or two regions are merged
00090     da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
00091 
00092     ++xs.x;
00093     ++xd.x;
00094     for(x = 1; x != w; ++x, ++xs.x, ++xd.x)
00095     {
00096         if((sa(xs) & Neighborhood::directionBit(Neighborhood::West)) ||
00097            (sa(xs, Neighborhood::west()) & Neighborhood::directionBit(Neighborhood::East)))
00098         {
00099             da.set(da(xd, Neighborhood::west()), xd);
00100         }
00101         else
00102         {
00103             da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
00104         }
00105     }
00106 
00107     ++ys.y;
00108     ++yd.y;
00109     for(y = 1; y != h; ++y, ++ys.y, ++yd.y)
00110     {
00111         xs = ys;
00112         xd = yd;
00113 
00114         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00115         {
00116             NeighborOffsetCirculator<Neighborhood> nc(x == w-1
00117                                                         ? ncstartBorder
00118                                                         : ncstart);
00119             NeighborOffsetCirculator<Neighborhood> nce(x == 0
00120                                                          ? ncendBorder
00121                                                          : ncend);
00122             LabelType currentLabel = labels.nextFreeLabel();
00123             for(; nc != nce; ++nc)
00124             {
00125                 if((sa(xs) & nc.directionBit()) || (sa(xs, *nc) & nc.oppositeDirectionBit()))
00126                 {
00127                     currentLabel = labels.makeUnion(da(xd,*nc), currentLabel);
00128                 }
00129             }
00130             da.set(labels.finalizeLabel(currentLabel), xd);
00131         }
00132     }
00133 
00134     unsigned int count = labels.makeContiguous();
00135 
00136     // pass 2: assign one label to each region (tree)
00137     // so that labels form a consecutive sequence 1, 2, ...
00138     yd = upperleftd;
00139     for(y=0; y != h; ++y, ++yd.y)
00140     {
00141         DestIterator xd(yd);
00142         for(x = 0; x != w; ++x, ++xd.x)
00143         {
00144             da.set(labels[da(xd)], xd);
00145         }
00146     }
00147     return count;
00148 }
00149 
00150 template <class SrcIterator, class SrcAccessor,
00151           class DestIterator, class DestAccessor,
00152           class Neighborhood>
00153 unsigned int watershedLabeling(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00154                                pair<DestIterator, DestAccessor> dest,
00155                                Neighborhood neighborhood)
00156 {
00157     return watershedLabeling(src.first, src.second, src.third,
00158                              dest.first, dest.second, neighborhood);
00159 }
00160 
00161 
00162 template <class SrcIterator, class SrcAccessor,
00163           class DestIterator, class DestAccessor>
00164 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00165                       DestIterator upperleftd, DestAccessor da,
00166                       FourNeighborCode)
00167 {
00168     int w = lowerrights.x - upperlefts.x;
00169     int h = lowerrights.y - upperlefts.y;
00170     int x,y;
00171 
00172     SrcIterator ys(upperlefts);
00173     SrcIterator xs(ys);
00174 
00175     DestIterator yd = upperleftd;
00176 
00177     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00178     {
00179         xs = ys;
00180         DestIterator xd = yd;
00181 
00182         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00183         {
00184             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00185             typename SrcAccessor::value_type v = sa(xs);
00186             // the following choice causes minima to point
00187             // to their lowest neighbor -- would this be better???
00188             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00189             int o = 0; // means center is minimum
00190             if(atBorder == NotAtBorder)
00191             {
00192                 NeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs), cend(c);
00193                 do {
00194                     if(sa(c) <= v)
00195                     {
00196                         v = sa(c);
00197                         o = c.directionBit();
00198                     }
00199                 }
00200                 while(++c != cend);
00201             }
00202             else
00203             {
00204                 RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs, atBorder), cend(c);
00205                 do {
00206                     if(sa(c) <= v)
00207                     {
00208                         v = sa(c);
00209                         o = c.directionBit();
00210                     }
00211                 }
00212                 while(++c != cend);
00213             }
00214             da.set(o, xd);
00215         }
00216     }
00217 }
00218 
00219 template <class SrcIterator, class SrcAccessor,
00220           class DestIterator, class DestAccessor>
00221 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00222                       DestIterator upperleftd, DestAccessor da,
00223                       EightNeighborCode)
00224 {
00225     int w = lowerrights.x - upperlefts.x;
00226     int h = lowerrights.y - upperlefts.y;
00227     int x,y;
00228 
00229     SrcIterator ys(upperlefts);
00230     SrcIterator xs(ys);
00231 
00232     DestIterator yd = upperleftd;
00233 
00234     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00235     {
00236         xs = ys;
00237         DestIterator xd = yd;
00238 
00239         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00240         {
00241             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00242             typename SrcAccessor::value_type v = sa(xs);
00243             // the following choice causes minima to point
00244             // to their lowest neighbor -- would this be better???
00245             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00246             int o = 0; // means center is minimum
00247             if(atBorder == NotAtBorder)
00248             {
00249                 // handle diagonal and principal neighbors separately
00250                 // so that principal neighbors are preferred when there are
00251                 // candidates with equal strength
00252                 NeighborhoodCirculator<SrcIterator, EightNeighborCode>
00253                                       c(xs, EightNeighborCode::NorthEast);
00254                 for(int i = 0; i < 4; ++i, c += 2)
00255                 {
00256                     if(sa(c) <= v)
00257                     {
00258                         v = sa(c);
00259                         o = c.directionBit();
00260                     }
00261                 }
00262                 --c;
00263                 for(int i = 0; i < 4; ++i, c += 2)
00264                 {
00265                     if(sa(c) <= v)
00266                     {
00267                         v = sa(c);
00268                         o = c.directionBit();
00269                     }
00270                 }
00271             }
00272             else
00273             {
00274                 RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode>
00275                              c(xs, atBorder), cend(c);
00276                 do
00277                 {
00278                     if(!c.isDiagonal())
00279                         continue;
00280                     if(sa(c) <= v)
00281                     {
00282                         v = sa(c);
00283                         o = c.directionBit();
00284                     }
00285                 }
00286                 while(++c != cend);
00287                 do
00288                 {
00289                     if(c.isDiagonal())
00290                         continue;
00291                     if(sa(c) <= v)
00292                     {
00293                         v = sa(c);
00294                         o = c.directionBit();
00295                     }
00296                 }
00297                 while(++c != cend);
00298             }
00299             da.set(o, xd);
00300         }
00301     }
00302 }
00303 
00304 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00305     Region growing, watersheds, and voronoi tesselation
00306 */
00307 //@{
00308 
00309 /********************************************************/
00310 /*                                                      */
00311 /*                      watersheds                      */
00312 /*                                                      */
00313 /********************************************************/
00314 
00315 /** \brief Region Segmentation by means of the watershed algorithm.
00316 
00317     This function implements the union-find version of the watershed algorithms
00318     as described in
00319 
00320     J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
00321     and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000
00322 
00323     The source image is a boundary indicator such as the gradient magnitude
00324     of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
00325     are used as region seeds, and all other pixels are recursively assigned to the same
00326     region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
00327     \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
00328     are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
00329     The function uses accessors.
00330 
00331     Note that VIGRA provides an alternative implementaion of the watershed transform via
00332     \ref seededRegionGrowing(). It is slower, but handles plateaus better
00333     and allows to keep a one pixel wide boundary between regions.
00334 
00335     <b> Declarations:</b>
00336 
00337     pass arguments explicitly:
00338     \code
00339     namespace vigra {
00340         template <class SrcIterator, class SrcAccessor,
00341                   class DestIterator, class DestAccessor,
00342                   class Neighborhood = EightNeighborCode>
00343         unsigned int
00344         watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00345                    DestIterator upperleftd, DestAccessor da,
00346                    Neighborhood neighborhood = EightNeighborCode())
00347     }
00348     \endcode
00349 
00350     use argument objects in conjunction with \ref ArgumentObjectFactories :
00351     \code
00352     namespace vigra {
00353         template <class SrcIterator, class SrcAccessor,
00354                   class DestIterator, class DestAccessor,
00355                   class Neighborhood = EightNeighborCode>
00356         unsigned int
00357         watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00358                    pair<DestIterator, DestAccessor> dest,
00359                    Neighborhood neighborhood = EightNeighborCode())
00360     }
00361     \endcode
00362 
00363     <b> Usage:</b>
00364 
00365     <b>\#include</b> <<a href="watersheds_8hxx-source.html">vigra/watersheds.hxx</a>><br>
00366     Namespace: vigra
00367 
00368     Example: watersheds of the gradient magnitude.
00369 
00370     \code
00371     vigra::BImage in(w,h);
00372     ... // read input data
00373 
00374     vigra::FImage gradx(x,y), grady(x,y), gradMag(x,y);
00375     gaussianGradient(srcImageRange(src), destImage(gradx), destImage(grady), 3.0);
00376     combineTwoImages(srcImageRange(gradx), srcImage(grady), destImage(gradMag),
00377                      vigra::MagnitudeFunctor<float>());
00378 
00379     // the pixel type of the destination image must be large enough to hold
00380     // numbers up to 'max_region_label' to prevent overflow
00381     vigra::IImage labeling(x,y);
00382     int max_region_label = watersheds(srcImageRange(gradMag), destImage(labeling));
00383 
00384     \endcode
00385 
00386     <b> Required Interface:</b>
00387 
00388     \code
00389     SrcImageIterator src_upperleft, src_lowerright;
00390     DestImageIterator dest_upperleft;
00391 
00392     SrcAccessor src_accessor;
00393     DestAccessor dest_accessor;
00394 
00395     // compare src values
00396     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
00397 
00398     // set result
00399     int label;
00400     dest_accessor.set(label, dest_upperleft);
00401     \endcode
00402 */
00403 doxygen_overloaded_function(template <...> unsigned int watersheds)
00404 
00405 template <class SrcIterator, class SrcAccessor,
00406           class DestIterator, class DestAccessor,
00407           class Neighborhood>
00408 unsigned int
00409 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00410            DestIterator upperleftd, DestAccessor da, Neighborhood neighborhood)
00411 {
00412     SImage orientationImage(lowerrights - upperlefts);
00413     SImage::traverser yo = orientationImage.upperLeft();
00414 
00415     prepareWatersheds(upperlefts, lowerrights, sa,
00416                      orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
00417     return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
00418                              upperleftd, da, neighborhood);
00419 }
00420 
00421 template <class SrcIterator, class SrcAccessor,
00422           class DestIterator, class DestAccessor>
00423 inline unsigned int
00424 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00425            DestIterator upperleftd, DestAccessor da)
00426 {
00427     return watersheds(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
00428 }
00429 
00430 template <class SrcIterator, class SrcAccessor,
00431           class DestIterator, class DestAccessor,
00432           class Neighborhood>
00433 inline unsigned int
00434 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00435            pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
00436 {
00437     return watersheds(src.first, src.second, src.third, dest.first, dest.second, neighborhood);
00438 }
00439 
00440 template <class SrcIterator, class SrcAccessor,
00441           class DestIterator, class DestAccessor>
00442 inline unsigned int
00443 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00444            pair<DestIterator, DestAccessor> dest)
00445 {
00446     return watersheds(src.first, src.second, src.third, dest.first, dest.second);
00447 }
00448 
00449 //@}
00450 
00451 } // namespace vigra
00452 
00453 #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.7.1 (Mon Apr 16 2012)