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

details vigra/watersheds.hxx VIGRA

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.5.0, Dec 07 2006 )                                    */
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 /*        koethe@informatik.uni-hamburg.de          or                  */
00012 /*        vigra@kogs1.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 template <class SrcIterator, class SrcAccessor,
00424           class DestIterator, class DestAccessor,
00425           class Neighborhood>
00426 unsigned int 
00427 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00428            DestIterator upperleftd, DestAccessor da, Neighborhood neighborhood)
00429 {
00430     SImage orientationImage(lowerrights - upperlefts);
00431     SImage::traverser yo = orientationImage.upperLeft();
00432 
00433     prepareWatersheds(upperlefts, lowerrights, sa, 
00434                      orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
00435     return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
00436                              upperleftd, da, neighborhood);
00437 }
00438 
00439 template <class SrcIterator, class SrcAccessor,
00440           class DestIterator, class DestAccessor>
00441 inline unsigned int 
00442 watersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
00443            DestIterator upperleftd, DestAccessor da)
00444 {
00445     return watersheds(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
00446 }
00447 
00448 template <class SrcIterator, class SrcAccessor,
00449           class DestIterator, class DestAccessor,
00450           class Neighborhood>
00451 inline unsigned int 
00452 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00453            pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
00454 {
00455     return watersheds(src.first, src.second, src.third, dest.first, dest.second, neighborhood);
00456 }
00457 
00458 template <class SrcIterator, class SrcAccessor,
00459           class DestIterator, class DestAccessor>
00460 inline unsigned int 
00461 watersheds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00462            pair<DestIterator, DestAccessor> dest)
00463 {
00464     return watersheds(src.first, src.second, src.third, dest.first, dest.second);
00465 }
00466 
00467 //@}
00468 
00469 } // namespace vigra
00470 
00471 #endif // VIGRA_WATERSHEDS_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)