[ 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 (3 Dec 2010)