[ 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 "localminmax.hxx" 00044 #include "labelimage.hxx" 00045 #include "seededregiongrowing.hxx" 00046 #include "functorexpression.hxx" 00047 #include "union_find.hxx" 00048 00049 namespace vigra { 00050 00051 template <class SrcIterator, class SrcAccessor, 00052 class DestIterator, class DestAccessor, 00053 class Neighborhood> 00054 unsigned int watershedLabeling(SrcIterator upperlefts, 00055 SrcIterator lowerrights, SrcAccessor sa, 00056 DestIterator upperleftd, DestAccessor da, 00057 Neighborhood) 00058 { 00059 typedef typename DestAccessor::value_type LabelType; 00060 00061 int w = lowerrights.x - upperlefts.x; 00062 int h = lowerrights.y - upperlefts.y; 00063 int x,y; 00064 00065 SrcIterator ys(upperlefts); 00066 SrcIterator xs(ys); 00067 DestIterator yd(upperleftd); 00068 DestIterator xd(yd); 00069 00070 // temporary image to store region labels 00071 detail::UnionFindArray<LabelType> labels; 00072 00073 // initialize the neighborhood circulators 00074 NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst); 00075 NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North); 00076 NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast); 00077 ++ncend; 00078 NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North); 00079 ++ncendBorder; 00080 00081 // pass 1: scan image from upper left to lower right 00082 // to find connected components 00083 00084 // Each component will be represented by a tree of pixels. Each 00085 // pixel contains the scan order address of its parent in the 00086 // tree. In order for pass 2 to work correctly, the parent must 00087 // always have a smaller scan order address than the child. 00088 // Therefore, we can merge trees only at their roots, because the 00089 // root of the combined tree must have the smallest scan order 00090 // address among all the tree's pixels/ nodes. The root of each 00091 // tree is distinguished by pointing to itself (it contains its 00092 // own scan order address). This condition is enforced whenever a 00093 // new region is found or two regions are merged 00094 da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd); 00095 00096 ++xs.x; 00097 ++xd.x; 00098 for(x = 1; x != w; ++x, ++xs.x, ++xd.x) 00099 { 00100 if((sa(xs) & Neighborhood::directionBit(Neighborhood::West)) || 00101 (sa(xs, Neighborhood::west()) & Neighborhood::directionBit(Neighborhood::East))) 00102 { 00103 da.set(da(xd, Neighborhood::west()), xd); 00104 } 00105 else 00106 { 00107 da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd); 00108 } 00109 } 00110 00111 ++ys.y; 00112 ++yd.y; 00113 for(y = 1; y != h; ++y, ++ys.y, ++yd.y) 00114 { 00115 xs = ys; 00116 xd = yd; 00117 00118 for(x = 0; x != w; ++x, ++xs.x, ++xd.x) 00119 { 00120 NeighborOffsetCirculator<Neighborhood> nc(x == w-1 00121 ? ncstartBorder 00122 : ncstart); 00123 NeighborOffsetCirculator<Neighborhood> nce(x == 0 00124 ? ncendBorder 00125 : ncend); 00126 LabelType currentLabel = labels.nextFreeLabel(); 00127 for(; nc != nce; ++nc) 00128 { 00129 if((sa(xs) & nc.directionBit()) || (sa(xs, *nc) & nc.oppositeDirectionBit())) 00130 { 00131 currentLabel = labels.makeUnion(da(xd,*nc), currentLabel); 00132 } 00133 } 00134 da.set(labels.finalizeLabel(currentLabel), xd); 00135 } 00136 } 00137 00138 unsigned int count = labels.makeContiguous(); 00139 00140 // pass 2: assign one label to each region (tree) 00141 // so that labels form a consecutive sequence 1, 2, ... 00142 yd = upperleftd; 00143 for(y=0; y != h; ++y, ++yd.y) 00144 { 00145 DestIterator xd(yd); 00146 for(x = 0; x != w; ++x, ++xd.x) 00147 { 00148 da.set(labels[da(xd)], xd); 00149 } 00150 } 00151 return count; 00152 } 00153 00154 template <class SrcIterator, class SrcAccessor, 00155 class DestIterator, class DestAccessor, 00156 class Neighborhood> 00157 unsigned int watershedLabeling(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00158 pair<DestIterator, DestAccessor> dest, 00159 Neighborhood neighborhood) 00160 { 00161 return watershedLabeling(src.first, src.second, src.third, 00162 dest.first, dest.second, neighborhood); 00163 } 00164 00165 00166 template <class SrcIterator, class SrcAccessor, 00167 class DestIterator, class DestAccessor> 00168 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00169 DestIterator upperleftd, DestAccessor da, 00170 FourNeighborCode) 00171 { 00172 int w = lowerrights.x - upperlefts.x; 00173 int h = lowerrights.y - upperlefts.y; 00174 int x,y; 00175 00176 SrcIterator ys(upperlefts); 00177 SrcIterator xs(ys); 00178 00179 DestIterator yd = upperleftd; 00180 00181 for(y = 0; y != h; ++y, ++ys.y, ++yd.y) 00182 { 00183 xs = ys; 00184 DestIterator xd = yd; 00185 00186 for(x = 0; x != w; ++x, ++xs.x, ++xd.x) 00187 { 00188 AtImageBorder atBorder = isAtImageBorder(x,y,w,h); 00189 typename SrcAccessor::value_type v = sa(xs); 00190 // the following choice causes minima to point 00191 // to their lowest neighbor -- would this be better??? 00192 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max(); 00193 int o = 0; // means center is minimum 00194 if(atBorder == NotAtBorder) 00195 { 00196 NeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs), cend(c); 00197 do { 00198 if(sa(c) <= v) 00199 { 00200 v = sa(c); 00201 o = c.directionBit(); 00202 } 00203 } 00204 while(++c != cend); 00205 } 00206 else 00207 { 00208 RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs, atBorder), cend(c); 00209 do { 00210 if(sa(c) <= v) 00211 { 00212 v = sa(c); 00213 o = c.directionBit(); 00214 } 00215 } 00216 while(++c != cend); 00217 } 00218 da.set(o, xd); 00219 } 00220 } 00221 } 00222 00223 template <class SrcIterator, class SrcAccessor, 00224 class DestIterator, class DestAccessor> 00225 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00226 DestIterator upperleftd, DestAccessor da, 00227 EightNeighborCode) 00228 { 00229 int w = lowerrights.x - upperlefts.x; 00230 int h = lowerrights.y - upperlefts.y; 00231 int x,y; 00232 00233 SrcIterator ys(upperlefts); 00234 SrcIterator xs(ys); 00235 00236 DestIterator yd = upperleftd; 00237 00238 for(y = 0; y != h; ++y, ++ys.y, ++yd.y) 00239 { 00240 xs = ys; 00241 DestIterator xd = yd; 00242 00243 for(x = 0; x != w; ++x, ++xs.x, ++xd.x) 00244 { 00245 AtImageBorder atBorder = isAtImageBorder(x,y,w,h); 00246 typename SrcAccessor::value_type v = sa(xs); 00247 // the following choice causes minima to point 00248 // to their lowest neighbor -- would this be better??? 00249 // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max(); 00250 int o = 0; // means center is minimum 00251 if(atBorder == NotAtBorder) 00252 { 00253 // handle diagonal and principal neighbors separately 00254 // so that principal neighbors are preferred when there are 00255 // candidates with equal strength 00256 NeighborhoodCirculator<SrcIterator, EightNeighborCode> 00257 c(xs, EightNeighborCode::NorthEast); 00258 for(int i = 0; i < 4; ++i, c += 2) 00259 { 00260 if(sa(c) <= v) 00261 { 00262 v = sa(c); 00263 o = c.directionBit(); 00264 } 00265 } 00266 --c; 00267 for(int i = 0; i < 4; ++i, c += 2) 00268 { 00269 if(sa(c) <= v) 00270 { 00271 v = sa(c); 00272 o = c.directionBit(); 00273 } 00274 } 00275 } 00276 else 00277 { 00278 RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode> 00279 c(xs, atBorder), cend(c); 00280 do 00281 { 00282 if(!c.isDiagonal()) 00283 continue; 00284 if(sa(c) <= v) 00285 { 00286 v = sa(c); 00287 o = c.directionBit(); 00288 } 00289 } 00290 while(++c != cend); 00291 do 00292 { 00293 if(c.isDiagonal()) 00294 continue; 00295 if(sa(c) <= v) 00296 { 00297 v = sa(c); 00298 o = c.directionBit(); 00299 } 00300 } 00301 while(++c != cend); 00302 } 00303 da.set(o, xd); 00304 } 00305 } 00306 } 00307 00308 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms 00309 Region growing, watersheds, and voronoi tesselation 00310 */ 00311 //@{ 00312 00313 /**\brief Options object for generateWatershedSeeds(). 00314 * 00315 <b> Usage:</b> 00316 00317 <b>\#include</b> <vigra/watersheds.hxx><br> 00318 Namespace: vigra 00319 00320 \code 00321 IImage seeds(boundary_indicator.size()); 00322 00323 // detect all minima in 'boundary_indicator' that are below gray level 22 00324 generateWatershedSeeds(srcImageRange(boundary_indicator), 00325 destImage(seeds), 00326 SeedOptions().minima().threshold(22.0)); 00327 \endcode 00328 */ 00329 class SeedOptions 00330 { 00331 public: 00332 enum DetectMinima { LevelSets, Minima, ExtendedMinima, Unspecified }; 00333 00334 double thresh; 00335 DetectMinima mini; 00336 00337 /**\brief Construct default options object. 00338 * 00339 Defaults are: detect minima without thresholding (i.e. all minima). 00340 */ 00341 SeedOptions() 00342 : thresh(NumericTraits<double>::max()), 00343 mini(Minima) 00344 {} 00345 00346 /** Generate seeds at minima. 00347 00348 Default: true 00349 */ 00350 SeedOptions & minima() 00351 { 00352 mini = Minima; 00353 return *this; 00354 } 00355 00356 /** Generate seeds at minima and minimal plateaus. 00357 00358 Default: false 00359 */ 00360 SeedOptions & extendedMinima() 00361 { 00362 mini = ExtendedMinima; 00363 return *this; 00364 } 00365 00366 /** Generate seeds as level sets. 00367 00368 Note that you must also set a threshold to define which level set is to be used.<br> 00369 Default: false 00370 */ 00371 SeedOptions & levelSets() 00372 { 00373 mini = LevelSets; 00374 return *this; 00375 } 00376 00377 /** Generate seeds as level sets at given threshold. 00378 00379 Equivalent to <tt>SeedOptions().levelSet().threshold(threshold)</tt><br> 00380 Default: false 00381 */ 00382 SeedOptions & levelSets(double threshold) 00383 { 00384 mini = LevelSets; 00385 thresh = threshold; 00386 return *this; 00387 } 00388 00389 /** Set threshold. 00390 00391 The threshold will be used by both the minima and level set variants 00392 of seed generation.<br> 00393 Default: no thresholding 00394 */ 00395 SeedOptions & threshold(double threshold) 00396 { 00397 thresh = threshold; 00398 return *this; 00399 } 00400 00401 // check whether the threshold has been set for the target type T 00402 template <class T> 00403 bool thresholdIsValid() const 00404 { 00405 return thresh < double(NumericTraits<T>::max()); 00406 } 00407 00408 // indicate that this option object is invalid (for internal use in watersheds) 00409 SeedOptions & unspecified() 00410 { 00411 mini = Unspecified; 00412 return *this; 00413 } 00414 }; 00415 00416 /** \brief Generate seeds for watershed computation and seeded region growing. 00417 00418 The source image is a boundary indicator such as the gradient magnitude 00419 or the trace of the \ref boundaryTensor(). Seeds are generally generated 00420 at locations where the boundaryness (i.e. the likelihood of the point being on the 00421 boundary) is very small. In particular, seeds can be placed by either 00422 looking for local minima (possibly including minimal plateaus) of the boundaryness, 00423 of by looking at level sets (i.e. regions where the boundaryness is below a threshold). 00424 Both methods can also be combined, so that only minima below a threshold are returned. 00425 The particular seeding strategy is specified by the <tt>options</tt> object 00426 (see \ref SeedOptions). 00427 00428 The pixel type of the input image must be <tt>LessThanComparable</tt>. 00429 The pixel type of the output image must be large enough to hold the labels for all seeds. 00430 (typically, you will use <tt>UInt32</tt>). The function will label seeds by consecutive integers 00431 (starting from 1) and returns the largest label it used. 00432 00433 Pass \ref vigra::EightNeighborCode or \ref vigra::FourNeighborCode to determine the 00434 neighborhood where pixel values are compared. 00435 00436 The function uses accessors. 00437 00438 <b> Declarations:</b> 00439 00440 pass arguments explicitly: 00441 \code 00442 namespace vigra { 00443 template <class SrcIterator, class SrcAccessor, 00444 class DestIterator, class DestAccessor, 00445 class Neighborhood = EightNeighborCode> 00446 unsigned int 00447 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00448 DestIterator upperleftd, DestAccessor da, 00449 Neighborhood neighborhood = EightNeighborCode(), 00450 SeedOptions const & options = SeedOptions()); 00451 } 00452 \endcode 00453 00454 use argument objects in conjunction with \ref ArgumentObjectFactories : 00455 \code 00456 namespace vigra { 00457 template <class SrcIterator, class SrcAccessor, 00458 class DestIterator, class DestAccessor, 00459 class Neighborhood = EightNeighborCode> 00460 unsigned int 00461 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00462 pair<DestIterator, DestAccessor> dest, 00463 Neighborhood neighborhood = EightNeighborCode(), 00464 SeedOptions const & options = SeedOptions()); 00465 } 00466 \endcode 00467 00468 <b> Usage:</b> 00469 00470 <b>\#include</b> <vigra/watersheds.hxx><br> 00471 Namespace: vigra 00472 00473 For detailed examples see watershedsRegionGrowing(). 00474 */ 00475 doxygen_overloaded_function(template <...> unsigned int generateWatershedSeeds) 00476 00477 template <class SrcIterator, class SrcAccessor, 00478 class DestIterator, class DestAccessor, 00479 class Neighborhood> 00480 unsigned int 00481 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00482 DestIterator upperleftd, DestAccessor da, 00483 Neighborhood neighborhood, 00484 SeedOptions const & options = SeedOptions()) 00485 { 00486 using namespace functor; 00487 typedef typename SrcAccessor::value_type SrcType; 00488 00489 vigra_precondition(options.mini != SeedOptions::LevelSets || 00490 options.thresholdIsValid<SrcType>(), 00491 "generateWatershedSeeds(): SeedOptions.levelSets() must be specified with threshold."); 00492 00493 Diff2D shape = lowerrights - upperlefts; 00494 BImage seeds(shape); 00495 00496 if(options.mini == SeedOptions::LevelSets) 00497 { 00498 transformImage(srcIterRange(upperlefts, lowerrights, sa), 00499 destImage(seeds), 00500 ifThenElse(Arg1() <= Param(options.thresh), Param(1), Param(0))); 00501 } 00502 else 00503 { 00504 localMinima(srcIterRange(upperlefts, lowerrights, sa), destImage(seeds), 00505 LocalMinmaxOptions().neighborhood(Neighborhood::DirectionCount) 00506 .markWith(1.0) 00507 .threshold(options.thresh) 00508 .allowAtBorder() 00509 .allowPlateaus(options.mini == SeedOptions::ExtendedMinima)); 00510 } 00511 00512 return labelImageWithBackground(srcImageRange(seeds), destIter(upperleftd, da), 00513 Neighborhood::DirectionCount == 8, 0); 00514 } 00515 00516 template <class SrcIterator, class SrcAccessor, 00517 class DestIterator, class DestAccessor> 00518 inline unsigned int 00519 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00520 DestIterator upperleftd, DestAccessor da, 00521 SeedOptions const & options = SeedOptions()) 00522 { 00523 return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da, 00524 EightNeighborCode(), options); 00525 } 00526 00527 template <class SrcIterator, class SrcAccessor, 00528 class DestIterator, class DestAccessor, 00529 class Neighborhood> 00530 inline unsigned int 00531 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00532 pair<DestIterator, DestAccessor> dest, 00533 Neighborhood neighborhood, 00534 SeedOptions const & options = SeedOptions()) 00535 { 00536 return generateWatershedSeeds(src.first, src.second, src.third, 00537 dest.first, dest.second, 00538 neighborhood, options); 00539 } 00540 00541 template <class SrcIterator, class SrcAccessor, 00542 class DestIterator, class DestAccessor> 00543 inline unsigned int 00544 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00545 pair<DestIterator, DestAccessor> dest, 00546 SeedOptions const & options = SeedOptions()) 00547 { 00548 return generateWatershedSeeds(src.first, src.second, src.third, 00549 dest.first, dest.second, 00550 EightNeighborCode(), options); 00551 } 00552 00553 /********************************************************/ 00554 /* */ 00555 /* watersheds */ 00556 /* */ 00557 /********************************************************/ 00558 00559 /** \brief Region segmentation by means of the union-find watershed algorithm. 00560 00561 This function implements the union-find version of the watershed algorithms 00562 as described in 00563 00564 J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms, 00565 and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000 00566 00567 The source image is a boundary indicator such as the gaussianGradientMagnitude() 00568 or the trace of the \ref boundaryTensor(). Local minima of the boundary indicator 00569 are used as region seeds, and all other pixels are recursively assigned to the same 00570 region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or 00571 \ref vigra::FourNeighborCode to determine the neighborhood where pixel values 00572 are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>. 00573 The function uses accessors. 00574 00575 Note that VIGRA provides an alternative implementation of the watershed transform via 00576 \ref watershedsRegionGrowing(). It is slower, but offers many more configuration options. 00577 00578 <b> Declarations:</b> 00579 00580 pass arguments explicitly: 00581 \code 00582 namespace vigra { 00583 template <class SrcIterator, class SrcAccessor, 00584 class DestIterator, class DestAccessor, 00585 class Neighborhood = EightNeighborCode> 00586 unsigned int 00587 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00588 DestIterator upperleftd, DestAccessor da, 00589 Neighborhood neighborhood = EightNeighborCode()) 00590 } 00591 \endcode 00592 00593 use argument objects in conjunction with \ref ArgumentObjectFactories : 00594 \code 00595 namespace vigra { 00596 template <class SrcIterator, class SrcAccessor, 00597 class DestIterator, class DestAccessor, 00598 class Neighborhood = EightNeighborCode> 00599 unsigned int 00600 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00601 pair<DestIterator, DestAccessor> dest, 00602 Neighborhood neighborhood = EightNeighborCode()) 00603 } 00604 \endcode 00605 00606 <b> Usage:</b> 00607 00608 <b>\#include</b> <vigra/watersheds.hxx><br> 00609 Namespace: vigra 00610 00611 Example: watersheds of the gradient magnitude. 00612 00613 \code 00614 vigra::BImage in(w,h); 00615 ... // read input data 00616 00617 // compute gradient magnitude as boundary indicator 00618 vigra::FImage gradMag(w, h); 00619 gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 3.0); 00620 00621 // the pixel type of the destination image must be large enough to hold 00622 // numbers up to 'max_region_label' to prevent overflow 00623 vigra::IImage labeling(w,h); 00624 int max_region_label = watershedsUnionFind(srcImageRange(gradMag), destImage(labeling)); 00625 00626 \endcode 00627 00628 <b> Required Interface:</b> 00629 00630 \code 00631 SrcIterator src_upperleft, src_lowerright; 00632 DestIterator dest_upperleft; 00633 00634 SrcAccessor src_accessor; 00635 DestAccessor dest_accessor; 00636 00637 // compare src values 00638 src_accessor(src_upperleft) <= src_accessor(src_upperleft) 00639 00640 // set result 00641 int label; 00642 dest_accessor.set(label, dest_upperleft); 00643 \endcode 00644 */ 00645 doxygen_overloaded_function(template <...> unsigned int watershedsUnionFind) 00646 00647 template <class SrcIterator, class SrcAccessor, 00648 class DestIterator, class DestAccessor, 00649 class Neighborhood> 00650 unsigned int 00651 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00652 DestIterator upperleftd, DestAccessor da, 00653 Neighborhood neighborhood) 00654 { 00655 SImage orientationImage(lowerrights - upperlefts); 00656 SImage::traverser yo = orientationImage.upperLeft(); 00657 00658 prepareWatersheds(upperlefts, lowerrights, sa, 00659 orientationImage.upperLeft(), orientationImage.accessor(), neighborhood); 00660 return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(), 00661 upperleftd, da, neighborhood); 00662 } 00663 00664 template <class SrcIterator, class SrcAccessor, 00665 class DestIterator, class DestAccessor> 00666 inline unsigned int 00667 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00668 DestIterator upperleftd, DestAccessor da) 00669 { 00670 return watershedsUnionFind(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode()); 00671 } 00672 00673 template <class SrcIterator, class SrcAccessor, 00674 class DestIterator, class DestAccessor, 00675 class Neighborhood> 00676 inline unsigned int 00677 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00678 pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood) 00679 { 00680 return watershedsUnionFind(src.first, src.second, src.third, 00681 dest.first, dest.second, neighborhood); 00682 } 00683 00684 template <class SrcIterator, class SrcAccessor, 00685 class DestIterator, class DestAccessor> 00686 inline unsigned int 00687 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00688 pair<DestIterator, DestAccessor> dest) 00689 { 00690 return watershedsUnionFind(src.first, src.second, src.third, 00691 dest.first, dest.second); 00692 } 00693 00694 /** \brief Options object for watershedsRegionGrowing(). 00695 00696 <b> Usage:</b> 00697 00698 see watershedsRegionGrowing() for detailed examples. 00699 */ 00700 class WatershedOptions 00701 { 00702 public: 00703 double max_cost, bias; 00704 SRGType terminate; 00705 unsigned int biased_label, bucket_count; 00706 SeedOptions seed_options; 00707 00708 /** \brief Create options object with default settings. 00709 00710 Defaults are: perform complete grow (all pixels are assigned to regions), 00711 use standard algorithm, assume that the destination image already contains 00712 region seeds. 00713 */ 00714 WatershedOptions() 00715 : max_cost(0.0), 00716 bias(1.0), 00717 terminate(CompleteGrow), 00718 biased_label(0), 00719 bucket_count(0), 00720 seed_options(SeedOptions().unspecified()) 00721 {} 00722 00723 /** \brief Perform complete grow. 00724 00725 That is, all pixels are assigned to regions, without explicit contours 00726 in between. 00727 00728 Default: true 00729 */ 00730 WatershedOptions & completeGrow() 00731 { 00732 terminate = SRGType(CompleteGrow | (terminate & StopAtThreshold)); 00733 return *this; 00734 } 00735 00736 /** \brief Keep one-pixel wide contour between regions. 00737 00738 Note that this option is unsupported by the turbo algorithm. 00739 00740 Default: false 00741 */ 00742 WatershedOptions & keepContours() 00743 { 00744 terminate = SRGType(KeepContours | (terminate & StopAtThreshold)); 00745 return *this; 00746 } 00747 00748 /** \brief Set \ref SRGType explicitly. 00749 00750 Default: CompleteGrow 00751 */ 00752 WatershedOptions & srgType(SRGType type) 00753 { 00754 terminate = type; 00755 return *this; 00756 } 00757 00758 /** \brief Stop region growing when the boundaryness exceeds the threshold. 00759 00760 This option may be combined with completeGrow() and keepContours(). 00761 00762 Default: no early stopping 00763 */ 00764 WatershedOptions & stopAtThreshold(double threshold) 00765 { 00766 terminate = SRGType(terminate | StopAtThreshold); 00767 max_cost = threshold; 00768 return *this; 00769 } 00770 00771 /** \brief Use a simpler, but faster region growing algorithm. 00772 00773 The algorithm internally uses a \ref BucketQueue to determine 00774 the processing order of the pixels. This is only useful, 00775 when the input boundary indicator image contains integers 00776 in the range <tt>[0, ..., bucket_count-1]</tt>. Since 00777 these boundary indicators are typically represented as 00778 UInt8 images, the default <tt>bucket_count</tt> is 256. 00779 00780 Default: don't use the turbo algorithm 00781 */ 00782 WatershedOptions & turboAlgorithm(unsigned int bucket_count = 256) 00783 { 00784 this->bucket_count = bucket_count; 00785 return *this; 00786 } 00787 00788 /** \brief Specify seed options. 00789 00790 In this case, watershedsRegionGrowing() assumes that the destination 00791 image does not yet contain seeds. It will therefore call 00792 generateWatershedSeeds() and pass on the seed options. 00793 00794 Default: don't compute seeds (i.e. assume that destination image already 00795 contains seeds). 00796 */ 00797 WatershedOptions & seedOptions(SeedOptions const & s) 00798 { 00799 seed_options = s; 00800 return *this; 00801 } 00802 00803 /** \brief Bias the cost of the specified region by the given factor. 00804 00805 In certain applications, one region (typically the background) should 00806 be preferred in region growing. This is most easily achieved 00807 by adjusting the assignment cost for that region as <tt>factor*cost</tt>, 00808 with a factor slightly below 1. 00809 00810 Default: don't bias any region. 00811 */ 00812 WatershedOptions & biasLabel(unsigned int label, double factor) 00813 { 00814 biased_label = label; 00815 bias = factor; 00816 return *this; 00817 } 00818 }; 00819 00820 namespace detail { 00821 00822 template <class CostType, class LabelType> 00823 class WatershedStatistics 00824 { 00825 public: 00826 00827 typedef SeedRgDirectValueFunctor<CostType> value_type; 00828 typedef value_type & reference; 00829 typedef value_type const & const_reference; 00830 00831 typedef CostType first_argument_type; 00832 typedef LabelType second_argument_type; 00833 typedef LabelType argument_type; 00834 00835 WatershedStatistics() 00836 {} 00837 00838 void resize(unsigned int) 00839 {} 00840 00841 void reset() 00842 {} 00843 00844 /** update regions statistics (do nothing in the watershed algorithm) 00845 */ 00846 template <class T1, class T2> 00847 void operator()(first_argument_type const &, second_argument_type const &) 00848 {} 00849 00850 /** ask for maximal index (label) allowed 00851 */ 00852 LabelType maxRegionLabel() const 00853 { return size() - 1; } 00854 00855 /** ask for array size (i.e. maxRegionLabel() + 1) 00856 */ 00857 LabelType size() const 00858 { return NumericTraits<LabelType>::max(); } 00859 00860 /** read the statistics functor for a region via its label 00861 */ 00862 const_reference operator[](argument_type label) const 00863 { return stats; } 00864 00865 /** access the statistics functor for a region via its label 00866 */ 00867 reference operator[](argument_type label) 00868 { return stats; } 00869 00870 value_type stats; 00871 }; 00872 00873 template <class Value> 00874 class SeedRgBiasedValueFunctor 00875 { 00876 public: 00877 double bias; 00878 00879 /* the functor's argument type 00880 */ 00881 typedef Value argument_type; 00882 00883 /* the functor's result type (unused, only necessary for 00884 use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics 00885 */ 00886 typedef Value result_type; 00887 00888 /* the return type of the cost() function 00889 */ 00890 typedef Value cost_type; 00891 00892 SeedRgBiasedValueFunctor(double b = 1.0) 00893 : bias(b) 00894 {} 00895 00896 /* Do nothing (since we need not update region statistics). 00897 */ 00898 void operator()(argument_type const &) const {} 00899 00900 /* Return scaled argument 00901 */ 00902 cost_type cost(argument_type const & v) const 00903 { 00904 return cost_type(bias*v); 00905 } 00906 }; 00907 00908 template <class CostType, class LabelType> 00909 class BiasedWatershedStatistics 00910 { 00911 public: 00912 00913 typedef SeedRgBiasedValueFunctor<CostType> value_type; 00914 typedef value_type & reference; 00915 typedef value_type const & const_reference; 00916 00917 typedef CostType first_argument_type; 00918 typedef LabelType second_argument_type; 00919 typedef LabelType argument_type; 00920 00921 BiasedWatershedStatistics(LabelType biasedLabel, double bias) 00922 : biased_label(biasedLabel), 00923 biased_stats(bias) 00924 {} 00925 00926 void resize(unsigned int) 00927 {} 00928 00929 void reset() 00930 {} 00931 00932 /** update regions statistics (do nothing in the watershed algorithm) 00933 */ 00934 template <class T1, class T2> 00935 void operator()(first_argument_type const &, second_argument_type const &) 00936 {} 00937 00938 /** ask for maximal index (label) allowed 00939 */ 00940 LabelType maxRegionLabel() const 00941 { return size() - 1; } 00942 00943 /** ask for array size (i.e. maxRegionLabel() + 1) 00944 */ 00945 LabelType size() const 00946 { return NumericTraits<LabelType>::max(); } 00947 00948 /** read the statistics functor for a region via its label 00949 */ 00950 const_reference operator[](argument_type label) const 00951 { 00952 return (label == biased_label) 00953 ? biased_stats 00954 : stats; 00955 } 00956 00957 /** access the statistics functor for a region via its label 00958 */ 00959 reference operator[](argument_type label) 00960 { 00961 return (label == biased_label) 00962 ? biased_stats 00963 : stats; 00964 } 00965 00966 LabelType biased_label; 00967 value_type stats, biased_stats; 00968 }; 00969 00970 } // namespace detail 00971 00972 /** \brief Region segmentation by means of a flooding-based watershed algorithm. 00973 00974 This function implements variants of the watershed algorithm 00975 described in 00976 00977 L. Vincent and P. Soille: "<em>Watersheds in digital spaces: An efficient algorithm 00978 based on immersion simulations</em>", IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991 00979 00980 The source image is a boundary indicator such as the gaussianGradientMagnitude() 00981 or the trace of the \ref boundaryTensor(), and the destination is a label image 00982 designating membership of each pixel in one of the regions. Plateaus in the boundary 00983 indicator (i.e. regions of constant gray value) are handled via a Euclidean distance 00984 transform by default. 00985 00986 By default, the destination image is assumed to hold seeds for a seeded watershed 00987 transform. Seeds may, for example, be created by means of generateWatershedSeeds(). 00988 Note that the seeds will be overridden with the final watershed segmentation. 00989 00990 Alternatively, you may provide \ref SeedOptions in order to instruct 00991 watershedsRegionGrowing() to generate its own seeds (it will call generateWatershedSeeds() 00992 internally). In that case, the destination image should be zero-initialized. 00993 00994 You can specify the neighborhood system to be used by passing \ref FourNeighborCode 00995 or \ref EightNeighborCode (default). 00996 00997 Further options to be specified via \ref WatershedOptions are: 00998 00999 <ul> 01000 <li> Whether to keep a 1-pixel-wide contour (with label 0) between regions or 01001 perform complete grow (i.e. all pixels are assigned to a region). 01002 <li> Whether to stop growing when the boundaryness exceeds a threshold (remaining 01003 pixels keep label 0). 01004 <li> Whether to use a faster, but less powerful algorithm ("turbo algorithm"). It 01005 is faster because it orders pixels by means of a \ref BucketQueue (therefore, 01006 the boundary indicator must contain integers in the range 01007 <tt>[0, ..., bucket_count-1]</tt>, where <tt>bucket_count</tt> is specified in 01008 the options object), it only supports complete growing (no contour between regions 01009 is possible), and it handles plateaus in a simplistic way. It also saves some 01010 memory because it allocates less temporary storage. 01011 <li> Whether one region (label) is to be preferred or discouraged by biasing its cost 01012 with a given factor (smaller than 1 for preference, larger than 1 for discouragement). 01013 </ul> 01014 01015 Note that VIGRA provides an alternative implementation of the watershed transform via 01016 \ref watershedsUnionFind(). 01017 01018 <b> Declarations:</b> 01019 01020 pass arguments explicitly: 01021 \code 01022 namespace vigra { 01023 template <class SrcIterator, class SrcAccessor, 01024 class DestIterator, class DestAccessor, 01025 class Neighborhood = EightNeighborCode> 01026 unsigned int 01027 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01028 DestIterator upperleftd, DestAccessor da, 01029 Neighborhood neighborhood = EightNeighborCode(), 01030 WatershedOptions const & options = WatershedOptions()); 01031 01032 template <class SrcIterator, class SrcAccessor, 01033 class DestIterator, class DestAccessor> 01034 unsigned int 01035 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01036 DestIterator upperleftd, DestAccessor da, 01037 WatershedOptions const & options = WatershedOptions()); 01038 } 01039 \endcode 01040 01041 use argument objects in conjunction with \ref ArgumentObjectFactories : 01042 \code 01043 namespace vigra { 01044 template <class SrcIterator, class SrcAccessor, 01045 class DestIterator, class DestAccessor, 01046 class Neighborhood = EightNeighborCode> 01047 unsigned int 01048 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01049 pair<DestIterator, DestAccessor> dest, 01050 Neighborhood neighborhood = EightNeighborCode(), 01051 WatershedOptions const & options = WatershedOptions()); 01052 01053 template <class SrcIterator, class SrcAccessor, 01054 class DestIterator, class DestAccessor> 01055 unsigned int 01056 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01057 pair<DestIterator, DestAccessor> dest, 01058 WatershedOptions const & options = WatershedOptions()); 01059 } 01060 \endcode 01061 01062 <b> Usage:</b> 01063 01064 <b>\#include</b> <vigra/watersheds.hxx><br> 01065 Namespace: vigra 01066 01067 Example: watersheds of the gradient magnitude. 01068 01069 \code 01070 vigra::BImage src(w, h); 01071 ... // read input data 01072 01073 // compute gradient magnitude at scale 1.0 as a boundary indicator 01074 vigra::FImage gradMag(w, h); 01075 gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0); 01076 01077 // example 1 01078 { 01079 // the pixel type of the destination image must be large enough to hold 01080 // numbers up to 'max_region_label' to prevent overflow 01081 vigra::IImage labeling(w, h); 01082 01083 // call watershed algorithm for 4-neighborhood, leave a 1-pixel boundary between regions, 01084 // and autogenerate seeds from all gradient minima where the magnitude is below 2.0 01085 unsigned int max_region_label = 01086 watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling), 01087 FourNeighborCode(), 01088 WatershedOptions().keepContours() 01089 .seedOptions(SeedOptions().minima().threshold(2.0))); 01090 } 01091 01092 // example 2 01093 { 01094 vigra::IImage labeling(w, h); 01095 01096 // compute seeds beforehand (use connected components of all pixels 01097 // where the gradient is below 4.0) 01098 unsigned int max_region_label = 01099 generateWatershedSeeds(srcImageRange(gradMag), destImage(labeling), 01100 SeedOptions().levelSets(4.0)); 01101 01102 // quantize the gradient image to 256 gray levels 01103 vigra::BImage gradMag256(w, h); 01104 vigra::FindMinMax<float> minmax; 01105 inspectImage(srcImageRange(gradMag), minmax); // find original range 01106 transformImage(srcImageRange(gradMag), destImage(gradMag256), 01107 linearRangeMapping(minmax, 0, 255)); 01108 01109 // call the turbo algorithm with 256 bins, using 8-neighborhood 01110 watershedsRegionGrowing(srcImageRange(gradMag256), destImage(labeling), 01111 WatershedOptions().turboAlgorithm(256)); 01112 } 01113 01114 // example 3 01115 { 01116 vigra::IImage labeling(w, h); 01117 01118 .. // get seeds from somewhere, e.g. an interactive labeling program, 01119 // make sure that label 1 corresponds to the background 01120 01121 // bias the watershed algorithm so that the background is preferred 01122 // by reducing the cost for label 1 to 90% 01123 watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling), 01124 WatershedOptions().biasLabel(1, 0.9)); 01125 } 01126 \endcode 01127 01128 <b> Required Interface:</b> 01129 01130 \code 01131 SrcIterator src_upperleft, src_lowerright; 01132 DestIterator dest_upperleft; 01133 01134 SrcAccessor src_accessor; 01135 DestAccessor dest_accessor; 01136 01137 // compare src values 01138 src_accessor(src_upperleft) <= src_accessor(src_upperleft) 01139 01140 // set result 01141 int label; 01142 dest_accessor.set(label, dest_upperleft); 01143 \endcode 01144 */ 01145 doxygen_overloaded_function(template <...> unsigned int watershedsRegionGrowing) 01146 01147 template <class SrcIterator, class SrcAccessor, 01148 class DestIterator, class DestAccessor, 01149 class Neighborhood> 01150 unsigned int 01151 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01152 DestIterator upperleftd, DestAccessor da, 01153 Neighborhood neighborhood, 01154 WatershedOptions const & options = WatershedOptions()) 01155 { 01156 typedef typename SrcAccessor::value_type ValueType; 01157 typedef typename DestAccessor::value_type LabelType; 01158 01159 unsigned int max_region_label = 0; 01160 01161 if(options.seed_options.mini != SeedOptions::Unspecified) 01162 { 01163 // we are supposed to compute seeds 01164 max_region_label = 01165 generateWatershedSeeds(srcIterRange(upperlefts, lowerrights, sa), 01166 destIter(upperleftd, da), 01167 neighborhood, options.seed_options); 01168 } 01169 01170 if(options.biased_label != 0) 01171 { 01172 // create a statistics functor for biased region growing 01173 detail::BiasedWatershedStatistics<ValueType, LabelType> 01174 regionstats(options.biased_label, options.bias); 01175 01176 // perform region growing, starting from the seeds computed above 01177 if(options.bucket_count == 0) 01178 { 01179 max_region_label = 01180 seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01181 srcIter(upperleftd, da), 01182 destIter(upperleftd, da), 01183 regionstats, options.terminate, neighborhood, options.max_cost); 01184 } 01185 else 01186 { 01187 max_region_label = 01188 fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01189 destIter(upperleftd, da), 01190 regionstats, options.terminate, 01191 neighborhood, options.max_cost, options.bucket_count); 01192 } 01193 } 01194 else 01195 { 01196 // create a statistics functor for region growing 01197 detail::WatershedStatistics<ValueType, LabelType> regionstats; 01198 01199 // perform region growing, starting from the seeds computed above 01200 if(options.bucket_count == 0) 01201 { 01202 max_region_label = 01203 seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01204 srcIter(upperleftd, da), 01205 destIter(upperleftd, da), 01206 regionstats, options.terminate, neighborhood, options.max_cost); 01207 } 01208 else 01209 { 01210 max_region_label = 01211 fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01212 destIter(upperleftd, da), 01213 regionstats, options.terminate, 01214 neighborhood, options.max_cost, options.bucket_count); 01215 } 01216 } 01217 01218 return max_region_label; 01219 } 01220 01221 template <class SrcIterator, class SrcAccessor, 01222 class DestIterator, class DestAccessor> 01223 inline unsigned int 01224 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01225 DestIterator upperleftd, DestAccessor da, 01226 WatershedOptions const & options = WatershedOptions()) 01227 { 01228 return watershedsRegionGrowing(upperlefts, lowerrights, sa, upperleftd, da, 01229 EightNeighborCode(), options); 01230 } 01231 01232 template <class SrcIterator, class SrcAccessor, 01233 class DestIterator, class DestAccessor, 01234 class Neighborhood> 01235 inline unsigned int 01236 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01237 pair<DestIterator, DestAccessor> dest, 01238 Neighborhood neighborhood, 01239 WatershedOptions const & options = WatershedOptions()) 01240 { 01241 return watershedsRegionGrowing(src.first, src.second, src.third, 01242 dest.first, dest.second, 01243 neighborhood, options); 01244 } 01245 01246 template <class SrcIterator, class SrcAccessor, 01247 class DestIterator, class DestAccessor> 01248 inline unsigned int 01249 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01250 pair<DestIterator, DestAccessor> dest, 01251 WatershedOptions const & options = WatershedOptions()) 01252 { 01253 return watershedsRegionGrowing(src.first, src.second, src.third, 01254 dest.first, dest.second, 01255 EightNeighborCode(), options); 01256 } 01257 01258 01259 //@} 01260 01261 } // namespace vigra 01262 01263 #endif // VIGRA_WATERSHEDS_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|