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