[ 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 LocalMinmaxOptions lm_options; 00505 lm_options.neighborhood(Neighborhood::DirectionCount) 00506 .markWith(1.0) 00507 .allowAtBorder() 00508 .allowPlateaus(options.mini == SeedOptions::ExtendedMinima); 00509 if(options.thresholdIsValid<SrcType>()) 00510 lm_options.threshold(options.thresh); 00511 00512 localMinima(srcIterRange(upperlefts, lowerrights, sa), destImage(seeds), 00513 lm_options); 00514 } 00515 00516 return labelImageWithBackground(srcImageRange(seeds), destIter(upperleftd, da), 00517 Neighborhood::DirectionCount == 8, 0); 00518 } 00519 00520 template <class SrcIterator, class SrcAccessor, 00521 class DestIterator, class DestAccessor> 00522 inline unsigned int 00523 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00524 DestIterator upperleftd, DestAccessor da, 00525 SeedOptions const & options = SeedOptions()) 00526 { 00527 return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da, 00528 EightNeighborCode(), options); 00529 } 00530 00531 template <class SrcIterator, class SrcAccessor, 00532 class DestIterator, class DestAccessor, 00533 class Neighborhood> 00534 inline unsigned int 00535 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00536 pair<DestIterator, DestAccessor> dest, 00537 Neighborhood neighborhood, 00538 SeedOptions const & options = SeedOptions()) 00539 { 00540 return generateWatershedSeeds(src.first, src.second, src.third, 00541 dest.first, dest.second, 00542 neighborhood, options); 00543 } 00544 00545 template <class SrcIterator, class SrcAccessor, 00546 class DestIterator, class DestAccessor> 00547 inline unsigned int 00548 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00549 pair<DestIterator, DestAccessor> dest, 00550 SeedOptions const & options = SeedOptions()) 00551 { 00552 return generateWatershedSeeds(src.first, src.second, src.third, 00553 dest.first, dest.second, 00554 EightNeighborCode(), options); 00555 } 00556 00557 /********************************************************/ 00558 /* */ 00559 /* watersheds */ 00560 /* */ 00561 /********************************************************/ 00562 00563 /** \brief Region segmentation by means of the union-find watershed algorithm. 00564 00565 This function implements the union-find version of the watershed algorithms 00566 as described in 00567 00568 J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms, 00569 and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000 00570 00571 The source image is a boundary indicator such as the gaussianGradientMagnitude() 00572 or the trace of the \ref boundaryTensor(). Local minima of the boundary indicator 00573 are used as region seeds, and all other pixels are recursively assigned to the same 00574 region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or 00575 \ref vigra::FourNeighborCode to determine the neighborhood where pixel values 00576 are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>. 00577 The function uses accessors. 00578 00579 Note that VIGRA provides an alternative implementation of the watershed transform via 00580 \ref watershedsRegionGrowing(). It is slower, but offers many more configuration options. 00581 00582 <b> Declarations:</b> 00583 00584 pass arguments explicitly: 00585 \code 00586 namespace vigra { 00587 template <class SrcIterator, class SrcAccessor, 00588 class DestIterator, class DestAccessor, 00589 class Neighborhood = EightNeighborCode> 00590 unsigned int 00591 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00592 DestIterator upperleftd, DestAccessor da, 00593 Neighborhood neighborhood = EightNeighborCode()) 00594 } 00595 \endcode 00596 00597 use argument objects in conjunction with \ref ArgumentObjectFactories : 00598 \code 00599 namespace vigra { 00600 template <class SrcIterator, class SrcAccessor, 00601 class DestIterator, class DestAccessor, 00602 class Neighborhood = EightNeighborCode> 00603 unsigned int 00604 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00605 pair<DestIterator, DestAccessor> dest, 00606 Neighborhood neighborhood = EightNeighborCode()) 00607 } 00608 \endcode 00609 00610 <b> Usage:</b> 00611 00612 <b>\#include</b> <vigra/watersheds.hxx><br> 00613 Namespace: vigra 00614 00615 Example: watersheds of the gradient magnitude. 00616 00617 \code 00618 vigra::BImage in(w,h); 00619 ... // read input data 00620 00621 // compute gradient magnitude as boundary indicator 00622 vigra::FImage gradMag(w, h); 00623 gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 3.0); 00624 00625 // the pixel type of the destination image must be large enough to hold 00626 // numbers up to 'max_region_label' to prevent overflow 00627 vigra::IImage labeling(w,h); 00628 int max_region_label = watershedsUnionFind(srcImageRange(gradMag), destImage(labeling)); 00629 00630 \endcode 00631 00632 <b> Required Interface:</b> 00633 00634 \code 00635 SrcIterator src_upperleft, src_lowerright; 00636 DestIterator dest_upperleft; 00637 00638 SrcAccessor src_accessor; 00639 DestAccessor dest_accessor; 00640 00641 // compare src values 00642 src_accessor(src_upperleft) <= src_accessor(src_upperleft) 00643 00644 // set result 00645 int label; 00646 dest_accessor.set(label, dest_upperleft); 00647 \endcode 00648 */ 00649 doxygen_overloaded_function(template <...> unsigned int watershedsUnionFind) 00650 00651 template <class SrcIterator, class SrcAccessor, 00652 class DestIterator, class DestAccessor, 00653 class Neighborhood> 00654 unsigned int 00655 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00656 DestIterator upperleftd, DestAccessor da, 00657 Neighborhood neighborhood) 00658 { 00659 SImage orientationImage(lowerrights - upperlefts); 00660 00661 prepareWatersheds(upperlefts, lowerrights, sa, 00662 orientationImage.upperLeft(), orientationImage.accessor(), neighborhood); 00663 return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(), 00664 upperleftd, da, neighborhood); 00665 } 00666 00667 template <class SrcIterator, class SrcAccessor, 00668 class DestIterator, class DestAccessor> 00669 inline unsigned int 00670 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 00671 DestIterator upperleftd, DestAccessor da) 00672 { 00673 return watershedsUnionFind(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode()); 00674 } 00675 00676 template <class SrcIterator, class SrcAccessor, 00677 class DestIterator, class DestAccessor, 00678 class Neighborhood> 00679 inline unsigned int 00680 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00681 pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood) 00682 { 00683 return watershedsUnionFind(src.first, src.second, src.third, 00684 dest.first, dest.second, neighborhood); 00685 } 00686 00687 template <class SrcIterator, class SrcAccessor, 00688 class DestIterator, class DestAccessor> 00689 inline unsigned int 00690 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src, 00691 pair<DestIterator, DestAccessor> dest) 00692 { 00693 return watershedsUnionFind(src.first, src.second, src.third, 00694 dest.first, dest.second); 00695 } 00696 00697 /** \brief Options object for watershedsRegionGrowing(). 00698 00699 <b> Usage:</b> 00700 00701 see watershedsRegionGrowing() for detailed examples. 00702 */ 00703 class WatershedOptions 00704 { 00705 public: 00706 double max_cost, bias; 00707 SRGType terminate; 00708 unsigned int biased_label, bucket_count; 00709 SeedOptions seed_options; 00710 00711 /** \brief Create options object with default settings. 00712 00713 Defaults are: perform complete grow (all pixels are assigned to regions), 00714 use standard algorithm, assume that the destination image already contains 00715 region seeds. 00716 */ 00717 WatershedOptions() 00718 : max_cost(0.0), 00719 bias(1.0), 00720 terminate(CompleteGrow), 00721 biased_label(0), 00722 bucket_count(0), 00723 seed_options(SeedOptions().unspecified()) 00724 {} 00725 00726 /** \brief Perform complete grow. 00727 00728 That is, all pixels are assigned to regions, without explicit contours 00729 in between. 00730 00731 Default: true 00732 */ 00733 WatershedOptions & completeGrow() 00734 { 00735 terminate = SRGType(CompleteGrow | (terminate & StopAtThreshold)); 00736 return *this; 00737 } 00738 00739 /** \brief Keep one-pixel wide contour between regions. 00740 00741 Note that this option is unsupported by the turbo algorithm. 00742 00743 Default: false 00744 */ 00745 WatershedOptions & keepContours() 00746 { 00747 terminate = SRGType(KeepContours | (terminate & StopAtThreshold)); 00748 return *this; 00749 } 00750 00751 /** \brief Set \ref SRGType explicitly. 00752 00753 Default: CompleteGrow 00754 */ 00755 WatershedOptions & srgType(SRGType type) 00756 { 00757 terminate = type; 00758 return *this; 00759 } 00760 00761 /** \brief Stop region growing when the boundaryness exceeds the threshold. 00762 00763 This option may be combined with completeGrow() and keepContours(). 00764 00765 Default: no early stopping 00766 */ 00767 WatershedOptions & stopAtThreshold(double threshold) 00768 { 00769 terminate = SRGType(terminate | StopAtThreshold); 00770 max_cost = threshold; 00771 return *this; 00772 } 00773 00774 /** \brief Use a simpler, but faster region growing algorithm. 00775 00776 The algorithm internally uses a \ref BucketQueue to determine 00777 the processing order of the pixels. This is only useful, 00778 when the input boundary indicator image contains integers 00779 in the range <tt>[0, ..., bucket_count-1]</tt>. Since 00780 these boundary indicators are typically represented as 00781 UInt8 images, the default <tt>bucket_count</tt> is 256. 00782 00783 Default: don't use the turbo algorithm 00784 */ 00785 WatershedOptions & turboAlgorithm(unsigned int bucket_count = 256) 00786 { 00787 this->bucket_count = bucket_count; 00788 return *this; 00789 } 00790 00791 /** \brief Specify seed options. 00792 00793 In this case, watershedsRegionGrowing() assumes that the destination 00794 image does not yet contain seeds. It will therefore call 00795 generateWatershedSeeds() and pass on the seed options. 00796 00797 Default: don't compute seeds (i.e. assume that destination image already 00798 contains seeds). 00799 */ 00800 WatershedOptions & seedOptions(SeedOptions const & s) 00801 { 00802 seed_options = s; 00803 return *this; 00804 } 00805 00806 /** \brief Bias the cost of the specified region by the given factor. 00807 00808 In certain applications, one region (typically the background) should 00809 be preferred in region growing. This is most easily achieved 00810 by adjusting the assignment cost for that region as <tt>factor*cost</tt>, 00811 with a factor slightly below 1. 00812 00813 Default: don't bias any region. 00814 */ 00815 WatershedOptions & biasLabel(unsigned int label, double factor) 00816 { 00817 biased_label = label; 00818 bias = factor; 00819 return *this; 00820 } 00821 }; 00822 00823 namespace detail { 00824 00825 template <class CostType, class LabelType> 00826 class WatershedStatistics 00827 { 00828 public: 00829 00830 typedef SeedRgDirectValueFunctor<CostType> value_type; 00831 typedef value_type & reference; 00832 typedef value_type const & const_reference; 00833 00834 typedef CostType first_argument_type; 00835 typedef LabelType second_argument_type; 00836 typedef LabelType argument_type; 00837 00838 WatershedStatistics() 00839 {} 00840 00841 void resize(unsigned int) 00842 {} 00843 00844 void reset() 00845 {} 00846 00847 /** update regions statistics (do nothing in the watershed algorithm) 00848 */ 00849 template <class T1, class T2> 00850 void operator()(first_argument_type const &, second_argument_type const &) 00851 {} 00852 00853 /** ask for maximal index (label) allowed 00854 */ 00855 LabelType maxRegionLabel() const 00856 { return size() - 1; } 00857 00858 /** ask for array size (i.e. maxRegionLabel() + 1) 00859 */ 00860 LabelType size() const 00861 { return NumericTraits<LabelType>::max(); } 00862 00863 /** read the statistics functor for a region via its label 00864 */ 00865 const_reference operator[](argument_type label) const 00866 { return stats; } 00867 00868 /** access the statistics functor for a region via its label 00869 */ 00870 reference operator[](argument_type label) 00871 { return stats; } 00872 00873 value_type stats; 00874 }; 00875 00876 template <class Value> 00877 class SeedRgBiasedValueFunctor 00878 { 00879 public: 00880 double bias; 00881 00882 /* the functor's argument type 00883 */ 00884 typedef Value argument_type; 00885 00886 /* the functor's result type (unused, only necessary for 00887 use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics 00888 */ 00889 typedef Value result_type; 00890 00891 /* the return type of the cost() function 00892 */ 00893 typedef Value cost_type; 00894 00895 SeedRgBiasedValueFunctor(double b = 1.0) 00896 : bias(b) 00897 {} 00898 00899 /* Do nothing (since we need not update region statistics). 00900 */ 00901 void operator()(argument_type const &) const {} 00902 00903 /* Return scaled argument 00904 */ 00905 cost_type cost(argument_type const & v) const 00906 { 00907 return cost_type(bias*v); 00908 } 00909 }; 00910 00911 template <class CostType, class LabelType> 00912 class BiasedWatershedStatistics 00913 { 00914 public: 00915 00916 typedef SeedRgBiasedValueFunctor<CostType> value_type; 00917 typedef value_type & reference; 00918 typedef value_type const & const_reference; 00919 00920 typedef CostType first_argument_type; 00921 typedef LabelType second_argument_type; 00922 typedef LabelType argument_type; 00923 00924 BiasedWatershedStatistics(LabelType biasedLabel, double bias) 00925 : biased_label(biasedLabel), 00926 biased_stats(bias) 00927 {} 00928 00929 void resize(unsigned int) 00930 {} 00931 00932 void reset() 00933 {} 00934 00935 /** update regions statistics (do nothing in the watershed algorithm) 00936 */ 00937 template <class T1, class T2> 00938 void operator()(first_argument_type const &, second_argument_type const &) 00939 {} 00940 00941 /** ask for maximal index (label) allowed 00942 */ 00943 LabelType maxRegionLabel() const 00944 { return size() - 1; } 00945 00946 /** ask for array size (i.e. maxRegionLabel() + 1) 00947 */ 00948 LabelType size() const 00949 { return NumericTraits<LabelType>::max(); } 00950 00951 /** read the statistics functor for a region via its label 00952 */ 00953 const_reference operator[](argument_type label) const 00954 { 00955 return (label == biased_label) 00956 ? biased_stats 00957 : stats; 00958 } 00959 00960 /** access the statistics functor for a region via its label 00961 */ 00962 reference operator[](argument_type label) 00963 { 00964 return (label == biased_label) 00965 ? biased_stats 00966 : stats; 00967 } 00968 00969 LabelType biased_label; 00970 value_type stats, biased_stats; 00971 }; 00972 00973 } // namespace detail 00974 00975 /** \brief Region segmentation by means of a flooding-based watershed algorithm. 00976 00977 This function implements variants of the watershed algorithm 00978 described in 00979 00980 L. Vincent and P. Soille: "<em>Watersheds in digital spaces: An efficient algorithm 00981 based on immersion simulations</em>", IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991 00982 00983 The source image is a boundary indicator such as the gaussianGradientMagnitude() 00984 or the trace of the \ref boundaryTensor(), and the destination is a label image 00985 designating membership of each pixel in one of the regions. Plateaus in the boundary 00986 indicator (i.e. regions of constant gray value) are handled via a Euclidean distance 00987 transform by default. 00988 00989 By default, the destination image is assumed to hold seeds for a seeded watershed 00990 transform. Seeds may, for example, be created by means of generateWatershedSeeds(). 00991 Note that the seeds will be overridden with the final watershed segmentation. 00992 00993 Alternatively, you may provide \ref SeedOptions in order to instruct 00994 watershedsRegionGrowing() to generate its own seeds (it will call generateWatershedSeeds() 00995 internally). In that case, the destination image should be zero-initialized. 00996 00997 You can specify the neighborhood system to be used by passing \ref FourNeighborCode 00998 or \ref EightNeighborCode (default). 00999 01000 Further options to be specified via \ref WatershedOptions are: 01001 01002 <ul> 01003 <li> Whether to keep a 1-pixel-wide contour (with label 0) between regions or 01004 perform complete grow (i.e. all pixels are assigned to a region). 01005 <li> Whether to stop growing when the boundaryness exceeds a threshold (remaining 01006 pixels keep label 0). 01007 <li> Whether to use a faster, but less powerful algorithm ("turbo algorithm"). It 01008 is faster because it orders pixels by means of a \ref BucketQueue (therefore, 01009 the boundary indicator must contain integers in the range 01010 <tt>[0, ..., bucket_count-1]</tt>, where <tt>bucket_count</tt> is specified in 01011 the options object), it only supports complete growing (no contour between regions 01012 is possible), and it handles plateaus in a simplistic way. It also saves some 01013 memory because it allocates less temporary storage. 01014 <li> Whether one region (label) is to be preferred or discouraged by biasing its cost 01015 with a given factor (smaller than 1 for preference, larger than 1 for discouragement). 01016 </ul> 01017 01018 Note that VIGRA provides an alternative implementation of the watershed transform via 01019 \ref watershedsUnionFind(). 01020 01021 <b> Declarations:</b> 01022 01023 pass arguments explicitly: 01024 \code 01025 namespace vigra { 01026 template <class SrcIterator, class SrcAccessor, 01027 class DestIterator, class DestAccessor, 01028 class Neighborhood = EightNeighborCode> 01029 unsigned int 01030 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01031 DestIterator upperleftd, DestAccessor da, 01032 Neighborhood neighborhood = EightNeighborCode(), 01033 WatershedOptions const & options = WatershedOptions()); 01034 01035 template <class SrcIterator, class SrcAccessor, 01036 class DestIterator, class DestAccessor> 01037 unsigned int 01038 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01039 DestIterator upperleftd, DestAccessor da, 01040 WatershedOptions const & options = WatershedOptions()); 01041 } 01042 \endcode 01043 01044 use argument objects in conjunction with \ref ArgumentObjectFactories : 01045 \code 01046 namespace vigra { 01047 template <class SrcIterator, class SrcAccessor, 01048 class DestIterator, class DestAccessor, 01049 class Neighborhood = EightNeighborCode> 01050 unsigned int 01051 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01052 pair<DestIterator, DestAccessor> dest, 01053 Neighborhood neighborhood = EightNeighborCode(), 01054 WatershedOptions const & options = WatershedOptions()); 01055 01056 template <class SrcIterator, class SrcAccessor, 01057 class DestIterator, class DestAccessor> 01058 unsigned int 01059 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01060 pair<DestIterator, DestAccessor> dest, 01061 WatershedOptions const & options = WatershedOptions()); 01062 } 01063 \endcode 01064 01065 <b> Usage:</b> 01066 01067 <b>\#include</b> <vigra/watersheds.hxx><br> 01068 Namespace: vigra 01069 01070 Example: watersheds of the gradient magnitude. 01071 01072 \code 01073 vigra::BImage src(w, h); 01074 ... // read input data 01075 01076 // compute gradient magnitude at scale 1.0 as a boundary indicator 01077 vigra::FImage gradMag(w, h); 01078 gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0); 01079 01080 // example 1 01081 { 01082 // the pixel type of the destination image must be large enough to hold 01083 // numbers up to 'max_region_label' to prevent overflow 01084 vigra::IImage labeling(w, h); 01085 01086 // call watershed algorithm for 4-neighborhood, leave a 1-pixel boundary between regions, 01087 // and autogenerate seeds from all gradient minima where the magnitude is below 2.0 01088 unsigned int max_region_label = 01089 watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling), 01090 FourNeighborCode(), 01091 WatershedOptions().keepContours() 01092 .seedOptions(SeedOptions().minima().threshold(2.0))); 01093 } 01094 01095 // example 2 01096 { 01097 vigra::IImage labeling(w, h); 01098 01099 // compute seeds beforehand (use connected components of all pixels 01100 // where the gradient is below 4.0) 01101 unsigned int max_region_label = 01102 generateWatershedSeeds(srcImageRange(gradMag), destImage(labeling), 01103 SeedOptions().levelSets(4.0)); 01104 01105 // quantize the gradient image to 256 gray levels 01106 vigra::BImage gradMag256(w, h); 01107 vigra::FindMinMax<float> minmax; 01108 inspectImage(srcImageRange(gradMag), minmax); // find original range 01109 transformImage(srcImageRange(gradMag), destImage(gradMag256), 01110 linearRangeMapping(minmax, 0, 255)); 01111 01112 // call the turbo algorithm with 256 bins, using 8-neighborhood 01113 watershedsRegionGrowing(srcImageRange(gradMag256), destImage(labeling), 01114 WatershedOptions().turboAlgorithm(256)); 01115 } 01116 01117 // example 3 01118 { 01119 vigra::IImage labeling(w, h); 01120 01121 .. // get seeds from somewhere, e.g. an interactive labeling program, 01122 // make sure that label 1 corresponds to the background 01123 01124 // bias the watershed algorithm so that the background is preferred 01125 // by reducing the cost for label 1 to 90% 01126 watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling), 01127 WatershedOptions().biasLabel(1, 0.9)); 01128 } 01129 \endcode 01130 01131 <b> Required Interface:</b> 01132 01133 \code 01134 SrcIterator src_upperleft, src_lowerright; 01135 DestIterator dest_upperleft; 01136 01137 SrcAccessor src_accessor; 01138 DestAccessor dest_accessor; 01139 01140 // compare src values 01141 src_accessor(src_upperleft) <= src_accessor(src_upperleft) 01142 01143 // set result 01144 int label; 01145 dest_accessor.set(label, dest_upperleft); 01146 \endcode 01147 */ 01148 doxygen_overloaded_function(template <...> unsigned int watershedsRegionGrowing) 01149 01150 template <class SrcIterator, class SrcAccessor, 01151 class DestIterator, class DestAccessor, 01152 class Neighborhood> 01153 unsigned int 01154 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01155 DestIterator upperleftd, DestAccessor da, 01156 Neighborhood neighborhood, 01157 WatershedOptions const & options = WatershedOptions()) 01158 { 01159 typedef typename SrcAccessor::value_type ValueType; 01160 typedef typename DestAccessor::value_type LabelType; 01161 01162 unsigned int max_region_label = 0; 01163 01164 if(options.seed_options.mini != SeedOptions::Unspecified) 01165 { 01166 // we are supposed to compute seeds 01167 max_region_label = 01168 generateWatershedSeeds(srcIterRange(upperlefts, lowerrights, sa), 01169 destIter(upperleftd, da), 01170 neighborhood, options.seed_options); 01171 } 01172 01173 if(options.biased_label != 0) 01174 { 01175 // create a statistics functor for biased region growing 01176 detail::BiasedWatershedStatistics<ValueType, LabelType> 01177 regionstats(options.biased_label, options.bias); 01178 01179 // perform region growing, starting from the seeds computed above 01180 if(options.bucket_count == 0) 01181 { 01182 max_region_label = 01183 seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01184 srcIter(upperleftd, da), 01185 destIter(upperleftd, da), 01186 regionstats, options.terminate, neighborhood, options.max_cost); 01187 } 01188 else 01189 { 01190 max_region_label = 01191 fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01192 destIter(upperleftd, da), 01193 regionstats, options.terminate, 01194 neighborhood, options.max_cost, options.bucket_count); 01195 } 01196 } 01197 else 01198 { 01199 // create a statistics functor for region growing 01200 detail::WatershedStatistics<ValueType, LabelType> regionstats; 01201 01202 // perform region growing, starting from the seeds computed above 01203 if(options.bucket_count == 0) 01204 { 01205 max_region_label = 01206 seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01207 srcIter(upperleftd, da), 01208 destIter(upperleftd, da), 01209 regionstats, options.terminate, neighborhood, options.max_cost); 01210 } 01211 else 01212 { 01213 max_region_label = 01214 fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa), 01215 destIter(upperleftd, da), 01216 regionstats, options.terminate, 01217 neighborhood, options.max_cost, options.bucket_count); 01218 } 01219 } 01220 01221 return max_region_label; 01222 } 01223 01224 template <class SrcIterator, class SrcAccessor, 01225 class DestIterator, class DestAccessor> 01226 inline unsigned int 01227 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa, 01228 DestIterator upperleftd, DestAccessor da, 01229 WatershedOptions const & options = WatershedOptions()) 01230 { 01231 return watershedsRegionGrowing(upperlefts, lowerrights, sa, upperleftd, da, 01232 EightNeighborCode(), options); 01233 } 01234 01235 template <class SrcIterator, class SrcAccessor, 01236 class DestIterator, class DestAccessor, 01237 class Neighborhood> 01238 inline unsigned int 01239 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01240 pair<DestIterator, DestAccessor> dest, 01241 Neighborhood neighborhood, 01242 WatershedOptions const & options = WatershedOptions()) 01243 { 01244 return watershedsRegionGrowing(src.first, src.second, src.third, 01245 dest.first, dest.second, 01246 neighborhood, options); 01247 } 01248 01249 template <class SrcIterator, class SrcAccessor, 01250 class DestIterator, class DestAccessor> 01251 inline unsigned int 01252 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01253 pair<DestIterator, DestAccessor> dest, 01254 WatershedOptions const & options = WatershedOptions()) 01255 { 01256 return watershedsRegionGrowing(src.first, src.second, src.third, 01257 dest.first, dest.second, 01258 EightNeighborCode(), options); 01259 } 01260 01261 01262 //@} 01263 01264 } // namespace vigra 01265 01266 #endif // VIGRA_WATERSHEDS_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|