[ 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 
00151 template <class SrcIterator, class SrcAccessor,
00152           class DestIterator, class DestAccessor>
00153 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00154                       DestIterator upperleftd, DestAccessor da,
00155                       FourNeighborCode)
00156 {
00157     int w = lowerrights.x - upperlefts.x;
00158     int h = lowerrights.y - upperlefts.y;
00159     int x,y;
00160 
00161     SrcIterator ys(upperlefts);
00162     SrcIterator xs(ys);
00163 
00164     DestIterator yd = upperleftd;
00165 
00166     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00167     {
00168         xs = ys;
00169         DestIterator xd = yd;
00170 
00171         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00172         {
00173             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00174             typename SrcAccessor::value_type v = sa(xs);
00175             // the following choice causes minima to point
00176             // to their lowest neighbor -- would this be better???
00177             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00178             int o = 0; // means center is minimum
00179             if(atBorder == NotAtBorder)
00180             {
00181                 NeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs), cend(c);
00182                 do {
00183                     if(sa(c) <= v)
00184                     {
00185                         v = sa(c);
00186                         o = c.directionBit();
00187                     }
00188                 }
00189                 while(++c != cend);
00190             }
00191             else
00192             {
00193                 RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode>  c(xs, atBorder), cend(c);
00194                 do {
00195                     if(sa(c) <= v)
00196                     {
00197                         v = sa(c);
00198                         o = c.directionBit();
00199                     }
00200                 }
00201                 while(++c != cend);
00202             }
00203             da.set(o, xd);
00204         }
00205     }
00206 }
00207 
00208 template <class SrcIterator, class SrcAccessor,
00209           class DestIterator, class DestAccessor>
00210 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00211                       DestIterator upperleftd, DestAccessor da,
00212                       EightNeighborCode)
00213 {
00214     int w = lowerrights.x - upperlefts.x;
00215     int h = lowerrights.y - upperlefts.y;
00216     int x,y;
00217 
00218     SrcIterator ys(upperlefts);
00219     SrcIterator xs(ys);
00220 
00221     DestIterator yd = upperleftd;
00222 
00223     for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
00224     {
00225         xs = ys;
00226         DestIterator xd = yd;
00227 
00228         for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
00229         {
00230             AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
00231             typename SrcAccessor::value_type v = sa(xs);
00232             // the following choice causes minima to point
00233             // to their lowest neighbor -- would this be better???
00234             // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
00235             int o = 0; // means center is minimum
00236             if(atBorder == NotAtBorder)
00237             {
00238                 // handle diagonal and principal neighbors separately
00239                 // so that principal neighbors are preferred when there are
00240                 // candidates with equal strength
00241                 NeighborhoodCirculator<SrcIterator, EightNeighborCode>
00242                                       c(xs, EightNeighborCode::NorthEast);
00243                 for(int i = 0; i < 4; ++i, c += 2)
00244                 {
00245                     if(sa(c) <= v)
00246                     {
00247                         v = sa(c);
00248                         o = c.directionBit();
00249                     }
00250                 }
00251                 --c;
00252                 for(int i = 0; i < 4; ++i, c += 2)
00253                 {
00254                     if(sa(c) <= v)
00255                     {
00256                         v = sa(c);
00257                         o = c.directionBit();
00258                     }
00259                 }
00260             }
00261             else
00262             {
00263                 RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode>
00264                              c(xs, atBorder), cend(c);
00265                 do
00266                 {
00267                     if(!c.isDiagonal())
00268                         continue;
00269                     if(sa(c) <= v)
00270                     {
00271                         v = sa(c);
00272                         o = c.directionBit();
00273                     }
00274                 }
00275                 while(++c != cend);
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             }
00288             da.set(o, xd);
00289         }
00290     }
00291 }
00292 
00293 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
00294     Region growing, watersheds, and voronoi tesselation
00295 */
00296 //@{
00297 
00298 /********************************************************/
00299 /*                                                      */
00300 /*                      watersheds                      */
00301 /*                                                      */
00302 /********************************************************/
00303 
00304 /** \brief Region Segmentation by means of the watershed algorithm.
00305 
00306     This function implements the union-find version of the watershed algorithms
00307     as described in
00308 
00309     J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
00310     and parallelization stretegies</em>", Fundamenta Informaticae, 41:187-228, 2000
00311 
00312     The source image is a boundary indicator such as the gradient magnitude
00313     of the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
00314     are used as region seeds, and all other pixels are recursively assigned to the same
00315     region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
00316     \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
00317     are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
00318     The function uses accessors.
00319 
00320     Note that VIGRA provides an alternative implementaion of the watershed transform via
00321     \ref seededRegionGrowing(). It is slower, but handles plateaus better
00322     and allows to keep a one pixel wide boundary between regions.
00323 
00324     <b> Declarations:</b>
00325 
00326     pass arguments explicitly:
00327     \code
00328     namespace vigra {
00329         template <class SrcIterator, class SrcAccessor,
00330                   class DestIterator, class DestAccessor,
00331                   class Neighborhood = EightNeighborCode>
00332         unsigned int
00333         watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00334                    DestIterator upperleftd, DestAccessor da,
00335                    Neighborhood neighborhood = EightNeighborCode())
00336     }
00337     \endcode
00338 
00339     use argument objects in conjunction with \ref ArgumentObjectFactories :
00340     \code
00341     namespace vigra {
00342         template <class SrcIterator, class SrcAccessor,
00343                   class DestIterator, class DestAccessor,
00344                   class Neighborhood = EightNeighborCode>
00345         unsigned int
00346         watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00347                    pair<DestIterator, DestAccessor> dest,
00348                    Neighborhood neighborhood = EightNeighborCode())
00349     }
00350     \endcode
00351 
00352     <b> Usage:</b>
00353 
00354     <b>\#include</b> <<a href="watersheds_8hxx-source.html">vigra/watersheds.hxx</a>><br>
00355     Namespace: vigra
00356 
00357     Example: watersheds of the gradient magnitude.
00358 
00359     \code
00360     vigra::BImage in(w,h);
00361     ... // read input data
00362 
00363     vigra::FImage gradx(x,y), grady(x,y), gradMag(x,y);
00364     gaussianGradient(srcImageRange(src), destImage(gradx), destImage(grady), 3.0);
00365     combineTwoImages(srcImageRange(gradx), srcImage(grady), destImage(gradMag),
00366                      vigra::MagnitudeFunctor<float>());
00367 
00368     // the pixel type of the destination image must be large enough to hold
00369     // numbers up to 'max_region_label' to prevent overflow
00370     vigra::IImage labeling(x,y);
00371     int max_region_label = watersheds(srcImageRange(gradMag), destImage(labeling));
00372 
00373     \endcode
00374 
00375     <b> Required Interface:</b>
00376 
00377     \code
00378     SrcImageIterator src_upperleft, src_lowerright;
00379     DestImageIterator dest_upperleft;
00380 
00381     SrcAccessor src_accessor;
00382     DestAccessor dest_accessor;
00383 
00384     // compare src values
00385     src_accessor(src_upperleft) <= src_accessor(src_upperleft)
00386 
00387     // set result
00388     int label;
00389     dest_accessor.set(label, dest_upperleft);
00390     \endcode
00391 */
00392 doxygen_overloaded_function(template <...> unsigned int watersheds)
00393 
00394 template <class SrcIterator, class SrcAccessor,
00395           class DestIterator, class DestAccessor,
00396           class Neighborhood>
00397 unsigned int
00398 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00399            DestIterator upperleftd, DestAccessor da, Neighborhood neighborhood)
00400 {
00401     SImage orientationImage(lowerrights - upperlefts);
00402     SImage::traverser yo = orientationImage.upperLeft();
00403 
00404     prepareWatersheds(upperlefts, lowerrights, sa,
00405                      orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
00406     return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
00407                              upperleftd, da, neighborhood);
00408 }
00409 
00410 template <class SrcIterator, class SrcAccessor,
00411           class DestIterator, class DestAccessor>
00412 inline unsigned int
00413 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00414            DestIterator upperleftd, DestAccessor da)
00415 {
00416     return watersheds(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
00417 }
00418 
00419 template <class SrcIterator, class SrcAccessor,
00420           class DestIterator, class DestAccessor,
00421           class Neighborhood>
00422 inline unsigned int
00423 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00424            pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
00425 {
00426     return watersheds(src.first, src.second, src.third, dest.first, dest.second, neighborhood);
00427 }
00428 
00429 template <class SrcIterator, class SrcAccessor,
00430           class DestIterator, class DestAccessor>
00431 inline unsigned int
00432 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00433            pair<DestIterator, DestAccessor> dest)
00434 {
00435     return watersheds(src.first, src.second, src.third, dest.first, dest.second);
00436 }
00437 
00438 //@}
00439 
00440 } // namespace vigra
00441 
00442 #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.0 (Thu Aug 25 2011)