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

watersheds.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2005 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_WATERSHEDS_HXX
37 #define VIGRA_WATERSHEDS_HXX
38 
39 #include <functional>
40 #include "mathutil.hxx"
41 #include "stdimage.hxx"
42 #include "pixelneighborhood.hxx"
43 #include "localminmax.hxx"
44 #include "labelimage.hxx"
45 #include "seededregiongrowing.hxx"
46 #include "functorexpression.hxx"
47 #include "union_find.hxx"
48 
49 namespace vigra {
50 
51 template <class SrcIterator, class SrcAccessor,
52  class DestIterator, class DestAccessor,
53  class Neighborhood>
54 unsigned int watershedLabeling(SrcIterator upperlefts,
55  SrcIterator lowerrights, SrcAccessor sa,
56  DestIterator upperleftd, DestAccessor da,
57  Neighborhood)
58 {
59  typedef typename DestAccessor::value_type LabelType;
60 
61  int w = lowerrights.x - upperlefts.x;
62  int h = lowerrights.y - upperlefts.y;
63  int x,y;
64 
65  SrcIterator ys(upperlefts);
66  SrcIterator xs(ys);
67  DestIterator yd(upperleftd);
68  DestIterator xd(yd);
69 
70  // temporary image to store region labels
71  detail::UnionFindArray<LabelType> labels;
72 
73  // initialize the neighborhood circulators
74  NeighborOffsetCirculator<Neighborhood> ncstart(Neighborhood::CausalFirst);
75  NeighborOffsetCirculator<Neighborhood> ncstartBorder(Neighborhood::North);
76  NeighborOffsetCirculator<Neighborhood> ncend(Neighborhood::CausalLast);
77  ++ncend;
78  NeighborOffsetCirculator<Neighborhood> ncendBorder(Neighborhood::North);
79  ++ncendBorder;
80 
81  // pass 1: scan image from upper left to lower right
82  // to find connected components
83 
84  // Each component will be represented by a tree of pixels. Each
85  // pixel contains the scan order address of its parent in the
86  // tree. In order for pass 2 to work correctly, the parent must
87  // always have a smaller scan order address than the child.
88  // Therefore, we can merge trees only at their roots, because the
89  // root of the combined tree must have the smallest scan order
90  // address among all the tree's pixels/ nodes. The root of each
91  // tree is distinguished by pointing to itself (it contains its
92  // own scan order address). This condition is enforced whenever a
93  // new region is found or two regions are merged
94  da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
95 
96  ++xs.x;
97  ++xd.x;
98  for(x = 1; x != w; ++x, ++xs.x, ++xd.x)
99  {
100  if((sa(xs) & Neighborhood::directionBit(Neighborhood::West)) ||
101  (sa(xs, Neighborhood::west()) & Neighborhood::directionBit(Neighborhood::East)))
102  {
103  da.set(da(xd, Neighborhood::west()), xd);
104  }
105  else
106  {
107  da.set(labels.finalizeLabel(labels.nextFreeLabel()), xd);
108  }
109  }
110 
111  ++ys.y;
112  ++yd.y;
113  for(y = 1; y != h; ++y, ++ys.y, ++yd.y)
114  {
115  xs = ys;
116  xd = yd;
117 
118  for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
119  {
120  NeighborOffsetCirculator<Neighborhood> nc(x == w-1
121  ? ncstartBorder
122  : ncstart);
123  NeighborOffsetCirculator<Neighborhood> nce(x == 0
124  ? ncendBorder
125  : ncend);
126  LabelType currentLabel = labels.nextFreeLabel();
127  for(; nc != nce; ++nc)
128  {
129  if((sa(xs) & nc.directionBit()) || (sa(xs, *nc) & nc.oppositeDirectionBit()))
130  {
131  currentLabel = labels.makeUnion(da(xd,*nc), currentLabel);
132  }
133  }
134  da.set(labels.finalizeLabel(currentLabel), xd);
135  }
136  }
137 
138  unsigned int count = labels.makeContiguous();
139 
140  // pass 2: assign one label to each region (tree)
141  // so that labels form a consecutive sequence 1, 2, ...
142  yd = upperleftd;
143  for(y=0; y != h; ++y, ++yd.y)
144  {
145  DestIterator xd(yd);
146  for(x = 0; x != w; ++x, ++xd.x)
147  {
148  da.set(labels[da(xd)], xd);
149  }
150  }
151  return count;
152 }
153 
154 template <class SrcIterator, class SrcAccessor,
155  class DestIterator, class DestAccessor,
156  class Neighborhood>
157 unsigned int watershedLabeling(triple<SrcIterator, SrcIterator, SrcAccessor> src,
158  pair<DestIterator, DestAccessor> dest,
159  Neighborhood neighborhood)
160 {
161  return watershedLabeling(src.first, src.second, src.third,
162  dest.first, dest.second, neighborhood);
163 }
164 
165 
166 template <class SrcIterator, class SrcAccessor,
167  class DestIterator, class DestAccessor>
168 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
169  DestIterator upperleftd, DestAccessor da,
171 {
172  int w = lowerrights.x - upperlefts.x;
173  int h = lowerrights.y - upperlefts.y;
174  int x,y;
175 
176  SrcIterator ys(upperlefts);
177  SrcIterator xs(ys);
178 
179  DestIterator yd = upperleftd;
180 
181  for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
182  {
183  xs = ys;
184  DestIterator xd = yd;
185 
186  for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
187  {
188  AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
189  typename SrcAccessor::value_type v = sa(xs);
190  // the following choice causes minima to point
191  // to their lowest neighbor -- would this be better???
192  // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
193  int o = 0; // means center is minimum
194  if(atBorder == NotAtBorder)
195  {
196  NeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs), cend(c);
197  do {
198  if(sa(c) <= v)
199  {
200  v = sa(c);
201  o = c.directionBit();
202  }
203  }
204  while(++c != cend);
205  }
206  else
207  {
208  RestrictedNeighborhoodCirculator<SrcIterator, FourNeighborCode> c(xs, atBorder), cend(c);
209  do {
210  if(sa(c) <= v)
211  {
212  v = sa(c);
213  o = c.directionBit();
214  }
215  }
216  while(++c != cend);
217  }
218  da.set(o, xd);
219  }
220  }
221 }
222 
223 template <class SrcIterator, class SrcAccessor,
224  class DestIterator, class DestAccessor>
225 void prepareWatersheds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
226  DestIterator upperleftd, DestAccessor da,
228 {
229  int w = lowerrights.x - upperlefts.x;
230  int h = lowerrights.y - upperlefts.y;
231  int x,y;
232 
233  SrcIterator ys(upperlefts);
234  SrcIterator xs(ys);
235 
236  DestIterator yd = upperleftd;
237 
238  for(y = 0; y != h; ++y, ++ys.y, ++yd.y)
239  {
240  xs = ys;
241  DestIterator xd = yd;
242 
243  for(x = 0; x != w; ++x, ++xs.x, ++xd.x)
244  {
245  AtImageBorder atBorder = isAtImageBorder(x,y,w,h);
246  typename SrcAccessor::value_type v = sa(xs);
247  // the following choice causes minima to point
248  // to their lowest neighbor -- would this be better???
249  // typename SrcAccessor::value_type v = NumericTraits<typename SrcAccessor::value_type>::max();
250  int o = 0; // means center is minimum
251  if(atBorder == NotAtBorder)
252  {
253  // handle diagonal and principal neighbors separately
254  // so that principal neighbors are preferred when there are
255  // candidates with equal strength
256  NeighborhoodCirculator<SrcIterator, EightNeighborCode>
258  for(int i = 0; i < 4; ++i, c += 2)
259  {
260  if(sa(c) <= v)
261  {
262  v = sa(c);
263  o = c.directionBit();
264  }
265  }
266  --c;
267  for(int i = 0; i < 4; ++i, c += 2)
268  {
269  if(sa(c) <= v)
270  {
271  v = sa(c);
272  o = c.directionBit();
273  }
274  }
275  }
276  else
277  {
278  RestrictedNeighborhoodCirculator<SrcIterator, EightNeighborCode>
279  c(xs, atBorder), cend(c);
280  do
281  {
282  if(!c.isDiagonal())
283  continue;
284  if(sa(c) <= v)
285  {
286  v = sa(c);
287  o = c.directionBit();
288  }
289  }
290  while(++c != cend);
291  do
292  {
293  if(c.isDiagonal())
294  continue;
295  if(sa(c) <= v)
296  {
297  v = sa(c);
298  o = c.directionBit();
299  }
300  }
301  while(++c != cend);
302  }
303  da.set(o, xd);
304  }
305  }
306 }
307 
308 /** \addtogroup SeededRegionGrowing Region Segmentation Algorithms
309  Region growing, watersheds, and voronoi tesselation
310 */
311 //@{
312 
313  /**\brief Options object for generateWatershedSeeds().
314  *
315  <b> Usage:</b>
316 
317  <b>\#include</b> <vigra/watersheds.hxx><br>
318  Namespace: vigra
319 
320  \code
321  IImage seeds(boundary_indicator.size());
322 
323  // detect all minima in 'boundary_indicator' that are below gray level 22
324  generateWatershedSeeds(srcImageRange(boundary_indicator),
325  destImage(seeds),
326  SeedOptions().minima().threshold(22.0));
327  \endcode
328  */
330 {
331 public:
332  enum DetectMinima { LevelSets, Minima, ExtendedMinima, Unspecified };
333 
334  double thresh;
335  DetectMinima mini;
336 
337  /**\brief Construct default options object.
338  *
339  Defaults are: detect minima without thresholding (i.e. all minima).
340  */
342  : thresh(NumericTraits<double>::max()),
343  mini(Minima)
344  {}
345 
346  /** Generate seeds at minima.
347 
348  Default: true
349  */
351  {
352  mini = Minima;
353  return *this;
354  }
355 
356  /** Generate seeds at minima and minimal plateaus.
357 
358  Default: false
359  */
361  {
362  mini = ExtendedMinima;
363  return *this;
364  }
365 
366  /** Generate seeds as level sets.
367 
368  Note that you must also set a threshold to define which level set is to be used.<br>
369  Default: false
370  */
372  {
373  mini = LevelSets;
374  return *this;
375  }
376 
377  /** Generate seeds as level sets at given threshold.
378 
379  Equivalent to <tt>SeedOptions().levelSet().threshold(threshold)</tt><br>
380  Default: false
381  */
383  {
384  mini = LevelSets;
385  thresh = threshold;
386  return *this;
387  }
388 
389  /** Set threshold.
390 
391  The threshold will be used by both the minima and level set variants
392  of seed generation.<br>
393  Default: no thresholding
394  */
396  {
397  thresh = threshold;
398  return *this;
399  }
400 
401  // check whether the threshold has been set for the target type T
402  template <class T>
403  bool thresholdIsValid() const
404  {
406  }
407 
408  // indicate that this option object is invalid (for internal use in watersheds)
409  SeedOptions & unspecified()
410  {
411  mini = Unspecified;
412  return *this;
413  }
414 };
415 
416 /** \brief Generate seeds for watershed computation and seeded region growing.
417 
418  The source image is a boundary indicator such as the gradient magnitude
419  or the trace of the \ref boundaryTensor(). Seeds are generally generated
420  at locations where the boundaryness (i.e. the likelihood of the point being on the
421  boundary) is very small. In particular, seeds can be placed by either
422  looking for local minima (possibly including minimal plateaus) of the boundaryness,
423  of by looking at level sets (i.e. regions where the boundaryness is below a threshold).
424  Both methods can also be combined, so that only minima below a threshold are returned.
425  The particular seeding strategy is specified by the <tt>options</tt> object
426  (see \ref SeedOptions).
427 
428  The pixel type of the input image must be <tt>LessThanComparable</tt>.
429  The pixel type of the output image must be large enough to hold the labels for all seeds.
430  (typically, you will use <tt>UInt32</tt>). The function will label seeds by consecutive integers
431  (starting from 1) and returns the largest label it used.
432 
433  Pass \ref vigra::EightNeighborCode or \ref vigra::FourNeighborCode to determine the
434  neighborhood where pixel values are compared.
435 
436  The function uses accessors.
437 
438  <b> Declarations:</b>
439 
440  pass arguments explicitly:
441  \code
442  namespace vigra {
443  template <class SrcIterator, class SrcAccessor,
444  class DestIterator, class DestAccessor,
445  class Neighborhood = EightNeighborCode>
446  unsigned int
447  generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
448  DestIterator upperleftd, DestAccessor da,
449  Neighborhood neighborhood = EightNeighborCode(),
450  SeedOptions const & options = SeedOptions());
451  }
452  \endcode
453 
454  use argument objects in conjunction with \ref ArgumentObjectFactories :
455  \code
456  namespace vigra {
457  template <class SrcIterator, class SrcAccessor,
458  class DestIterator, class DestAccessor,
459  class Neighborhood = EightNeighborCode>
460  unsigned int
461  generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
462  pair<DestIterator, DestAccessor> dest,
463  Neighborhood neighborhood = EightNeighborCode(),
464  SeedOptions const & options = SeedOptions());
465  }
466  \endcode
467 
468  <b> Usage:</b>
469 
470  <b>\#include</b> <vigra/watersheds.hxx><br>
471  Namespace: vigra
472 
473  For detailed examples see watershedsRegionGrowing().
474 */
475 doxygen_overloaded_function(template <...> unsigned int generateWatershedSeeds)
476 
477 template <class SrcIterator, class SrcAccessor,
478  class DestIterator, class DestAccessor,
479  class Neighborhood>
480 unsigned int
481 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
482  DestIterator upperleftd, DestAccessor da,
483  Neighborhood neighborhood,
484  SeedOptions const & options = SeedOptions())
485 {
486  using namespace functor;
487  typedef typename SrcAccessor::value_type SrcType;
488 
489  vigra_precondition(options.mini != SeedOptions::LevelSets ||
490  options.thresholdIsValid<SrcType>(),
491  "generateWatershedSeeds(): SeedOptions.levelSets() must be specified with threshold.");
492 
493  Diff2D shape = lowerrights - upperlefts;
494  BImage seeds(shape);
495 
496  if(options.mini == SeedOptions::LevelSets)
497  {
498  transformImage(srcIterRange(upperlefts, lowerrights, sa),
499  destImage(seeds),
500  ifThenElse(Arg1() <= Param(options.thresh), Param(1), Param(0)));
501  }
502  else
503  {
504  localMinima(srcIterRange(upperlefts, lowerrights, sa), destImage(seeds),
505  LocalMinmaxOptions().neighborhood(Neighborhood::DirectionCount)
506  .markWith(1.0)
507  .threshold(options.thresh)
508  .allowAtBorder()
509  .allowPlateaus(options.mini == SeedOptions::ExtendedMinima));
510  }
511 
512  return labelImageWithBackground(srcImageRange(seeds), destIter(upperleftd, da),
513  Neighborhood::DirectionCount == 8, 0);
514 }
515 
516 template <class SrcIterator, class SrcAccessor,
517  class DestIterator, class DestAccessor>
518 inline unsigned int
519 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
520  DestIterator upperleftd, DestAccessor da,
521  SeedOptions const & options = SeedOptions())
522 {
523  return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da,
524  EightNeighborCode(), options);
525 }
526 
527 template <class SrcIterator, class SrcAccessor,
528  class DestIterator, class DestAccessor,
529  class Neighborhood>
530 inline unsigned int
531 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
532  pair<DestIterator, DestAccessor> dest,
533  Neighborhood neighborhood,
534  SeedOptions const & options = SeedOptions())
535 {
536  return generateWatershedSeeds(src.first, src.second, src.third,
537  dest.first, dest.second,
538  neighborhood, options);
539 }
540 
541 template <class SrcIterator, class SrcAccessor,
542  class DestIterator, class DestAccessor>
543 inline unsigned int
544 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
545  pair<DestIterator, DestAccessor> dest,
546  SeedOptions const & options = SeedOptions())
547 {
548  return generateWatershedSeeds(src.first, src.second, src.third,
549  dest.first, dest.second,
550  EightNeighborCode(), options);
551 }
552 
553 /********************************************************/
554 /* */
555 /* watersheds */
556 /* */
557 /********************************************************/
558 
559 /** \brief Region segmentation by means of the union-find watershed algorithm.
560 
561  This function implements the union-find version of the watershed algorithms
562  as described in
563 
564  J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
565  and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000
566 
567  The source image is a boundary indicator such as the gaussianGradientMagnitude()
568  or the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
569  are used as region seeds, and all other pixels are recursively assigned to the same
570  region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
571  \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
572  are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
573  The function uses accessors.
574 
575  Note that VIGRA provides an alternative implementation of the watershed transform via
576  \ref watershedsRegionGrowing(). It is slower, but offers many more configuration options.
577 
578  <b> Declarations:</b>
579 
580  pass arguments explicitly:
581  \code
582  namespace vigra {
583  template <class SrcIterator, class SrcAccessor,
584  class DestIterator, class DestAccessor,
585  class Neighborhood = EightNeighborCode>
586  unsigned int
587  watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
588  DestIterator upperleftd, DestAccessor da,
589  Neighborhood neighborhood = EightNeighborCode())
590  }
591  \endcode
592 
593  use argument objects in conjunction with \ref ArgumentObjectFactories :
594  \code
595  namespace vigra {
596  template <class SrcIterator, class SrcAccessor,
597  class DestIterator, class DestAccessor,
598  class Neighborhood = EightNeighborCode>
599  unsigned int
600  watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
601  pair<DestIterator, DestAccessor> dest,
602  Neighborhood neighborhood = EightNeighborCode())
603  }
604  \endcode
605 
606  <b> Usage:</b>
607 
608  <b>\#include</b> <vigra/watersheds.hxx><br>
609  Namespace: vigra
610 
611  Example: watersheds of the gradient magnitude.
612 
613  \code
614  vigra::BImage in(w,h);
615  ... // read input data
616 
617  // compute gradient magnitude as boundary indicator
618  vigra::FImage gradMag(w, h);
619  gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 3.0);
620 
621  // the pixel type of the destination image must be large enough to hold
622  // numbers up to 'max_region_label' to prevent overflow
623  vigra::IImage labeling(w,h);
624  int max_region_label = watershedsUnionFind(srcImageRange(gradMag), destImage(labeling));
625 
626  \endcode
627 
628  <b> Required Interface:</b>
629 
630  \code
631  SrcIterator src_upperleft, src_lowerright;
632  DestIterator dest_upperleft;
633 
634  SrcAccessor src_accessor;
635  DestAccessor dest_accessor;
636 
637  // compare src values
638  src_accessor(src_upperleft) <= src_accessor(src_upperleft)
639 
640  // set result
641  int label;
642  dest_accessor.set(label, dest_upperleft);
643  \endcode
644 */
645 doxygen_overloaded_function(template <...> unsigned int watershedsUnionFind)
646 
647 template <class SrcIterator, class SrcAccessor,
648  class DestIterator, class DestAccessor,
649  class Neighborhood>
650 unsigned int
651 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
652  DestIterator upperleftd, DestAccessor da,
653  Neighborhood neighborhood)
654 {
655  SImage orientationImage(lowerrights - upperlefts);
656  SImage::traverser yo = orientationImage.upperLeft();
657 
658  prepareWatersheds(upperlefts, lowerrights, sa,
659  orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
660  return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
661  upperleftd, da, neighborhood);
662 }
663 
664 template <class SrcIterator, class SrcAccessor,
665  class DestIterator, class DestAccessor>
666 inline unsigned int
667 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
668  DestIterator upperleftd, DestAccessor da)
669 {
670  return watershedsUnionFind(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
671 }
672 
673 template <class SrcIterator, class SrcAccessor,
674  class DestIterator, class DestAccessor,
675  class Neighborhood>
676 inline unsigned int
677 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
678  pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
679 {
680  return watershedsUnionFind(src.first, src.second, src.third,
681  dest.first, dest.second, neighborhood);
682 }
683 
684 template <class SrcIterator, class SrcAccessor,
685  class DestIterator, class DestAccessor>
686 inline unsigned int
687 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
688  pair<DestIterator, DestAccessor> dest)
689 {
690  return watershedsUnionFind(src.first, src.second, src.third,
691  dest.first, dest.second);
692 }
693 
694 /** \brief Options object for watershedsRegionGrowing().
695 
696  <b> Usage:</b>
697 
698  see watershedsRegionGrowing() for detailed examples.
699 */
701 {
702  public:
703  double max_cost, bias;
704  SRGType terminate;
705  unsigned int biased_label, bucket_count;
706  SeedOptions seed_options;
707 
708  /** \brief Create options object with default settings.
709 
710  Defaults are: perform complete grow (all pixels are assigned to regions),
711  use standard algorithm, assume that the destination image already contains
712  region seeds.
713  */
715  : max_cost(0.0),
716  bias(1.0),
717  terminate(CompleteGrow),
718  biased_label(0),
719  bucket_count(0),
720  seed_options(SeedOptions().unspecified())
721  {}
722 
723  /** \brief Perform complete grow.
724 
725  That is, all pixels are assigned to regions, without explicit contours
726  in between.
727 
728  Default: true
729  */
731  {
732  terminate = SRGType(CompleteGrow | (terminate & StopAtThreshold));
733  return *this;
734  }
735 
736  /** \brief Keep one-pixel wide contour between regions.
737 
738  Note that this option is unsupported by the turbo algorithm.
739 
740  Default: false
741  */
743  {
744  terminate = SRGType(KeepContours | (terminate & StopAtThreshold));
745  return *this;
746  }
747 
748  /** \brief Set \ref SRGType explicitly.
749 
750  Default: CompleteGrow
751  */
753  {
754  terminate = type;
755  return *this;
756  }
757 
758  /** \brief Stop region growing when the boundaryness exceeds the threshold.
759 
760  This option may be combined with completeGrow() and keepContours().
761 
762  Default: no early stopping
763  */
764  WatershedOptions & stopAtThreshold(double threshold)
765  {
766  terminate = SRGType(terminate | StopAtThreshold);
767  max_cost = threshold;
768  return *this;
769  }
770 
771  /** \brief Use a simpler, but faster region growing algorithm.
772 
773  The algorithm internally uses a \ref BucketQueue to determine
774  the processing order of the pixels. This is only useful,
775  when the input boundary indicator image contains integers
776  in the range <tt>[0, ..., bucket_count-1]</tt>. Since
777  these boundary indicators are typically represented as
778  UInt8 images, the default <tt>bucket_count</tt> is 256.
779 
780  Default: don't use the turbo algorithm
781  */
782  WatershedOptions & turboAlgorithm(unsigned int bucket_count = 256)
783  {
784  this->bucket_count = bucket_count;
785  return *this;
786  }
787 
788  /** \brief Specify seed options.
789 
790  In this case, watershedsRegionGrowing() assumes that the destination
791  image does not yet contain seeds. It will therefore call
792  generateWatershedSeeds() and pass on the seed options.
793 
794  Default: don't compute seeds (i.e. assume that destination image already
795  contains seeds).
796  */
798  {
799  seed_options = s;
800  return *this;
801  }
802 
803  /** \brief Bias the cost of the specified region by the given factor.
804 
805  In certain applications, one region (typically the background) should
806  be preferred in region growing. This is most easily achieved
807  by adjusting the assignment cost for that region as <tt>factor*cost</tt>,
808  with a factor slightly below 1.
809 
810  Default: don't bias any region.
811  */
812  WatershedOptions & biasLabel(unsigned int label, double factor)
813  {
814  biased_label = label;
815  bias = factor;
816  return *this;
817  }
818 };
819 
820 namespace detail {
821 
822 template <class CostType, class LabelType>
823 class WatershedStatistics
824 {
825  public:
826 
827  typedef SeedRgDirectValueFunctor<CostType> value_type;
828  typedef value_type & reference;
829  typedef value_type const & const_reference;
830 
831  typedef CostType first_argument_type;
832  typedef LabelType second_argument_type;
833  typedef LabelType argument_type;
834 
835  WatershedStatistics()
836  {}
837 
838  void resize(unsigned int)
839  {}
840 
841  void reset()
842  {}
843 
844  /** update regions statistics (do nothing in the watershed algorithm)
845  */
846  template <class T1, class T2>
847  void operator()(first_argument_type const &, second_argument_type const &)
848  {}
849 
850  /** ask for maximal index (label) allowed
851  */
852  LabelType maxRegionLabel() const
853  { return size() - 1; }
854 
855  /** ask for array size (i.e. maxRegionLabel() + 1)
856  */
857  LabelType size() const
858  { return NumericTraits<LabelType>::max(); }
859 
860  /** read the statistics functor for a region via its label
861  */
862  const_reference operator[](argument_type label) const
863  { return stats; }
864 
865  /** access the statistics functor for a region via its label
866  */
867  reference operator[](argument_type label)
868  { return stats; }
869 
870  value_type stats;
871 };
872 
873 template <class Value>
874 class SeedRgBiasedValueFunctor
875 {
876  public:
877  double bias;
878 
879  /* the functor's argument type
880  */
881  typedef Value argument_type;
882 
883  /* the functor's result type (unused, only necessary for
884  use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
885  */
886  typedef Value result_type;
887 
888  /* the return type of the cost() function
889  */
890  typedef Value cost_type;
891 
892  SeedRgBiasedValueFunctor(double b = 1.0)
893  : bias(b)
894  {}
895 
896  /* Do nothing (since we need not update region statistics).
897  */
898  void operator()(argument_type const &) const {}
899 
900  /* Return scaled argument
901  */
902  cost_type cost(argument_type const & v) const
903  {
904  return cost_type(bias*v);
905  }
906 };
907 
908 template <class CostType, class LabelType>
909 class BiasedWatershedStatistics
910 {
911  public:
912 
913  typedef SeedRgBiasedValueFunctor<CostType> value_type;
914  typedef value_type & reference;
915  typedef value_type const & const_reference;
916 
917  typedef CostType first_argument_type;
918  typedef LabelType second_argument_type;
919  typedef LabelType argument_type;
920 
921  BiasedWatershedStatistics(LabelType biasedLabel, double bias)
922  : biased_label(biasedLabel),
923  biased_stats(bias)
924  {}
925 
926  void resize(unsigned int)
927  {}
928 
929  void reset()
930  {}
931 
932  /** update regions statistics (do nothing in the watershed algorithm)
933  */
934  template <class T1, class T2>
935  void operator()(first_argument_type const &, second_argument_type const &)
936  {}
937 
938  /** ask for maximal index (label) allowed
939  */
940  LabelType maxRegionLabel() const
941  { return size() - 1; }
942 
943  /** ask for array size (i.e. maxRegionLabel() + 1)
944  */
945  LabelType size() const
946  { return NumericTraits<LabelType>::max(); }
947 
948  /** read the statistics functor for a region via its label
949  */
950  const_reference operator[](argument_type label) const
951  {
952  return (label == biased_label)
953  ? biased_stats
954  : stats;
955  }
956 
957  /** access the statistics functor for a region via its label
958  */
959  reference operator[](argument_type label)
960  {
961  return (label == biased_label)
962  ? biased_stats
963  : stats;
964  }
965 
966  LabelType biased_label;
967  value_type stats, biased_stats;
968 };
969 
970 } // namespace detail
971 
972 /** \brief Region segmentation by means of a flooding-based watershed algorithm.
973 
974  This function implements variants of the watershed algorithm
975  described in
976 
977  L. Vincent and P. Soille: "<em>Watersheds in digital spaces: An efficient algorithm
978  based on immersion simulations</em>", IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991
979 
980  The source image is a boundary indicator such as the gaussianGradientMagnitude()
981  or the trace of the \ref boundaryTensor(), and the destination is a label image
982  designating membership of each pixel in one of the regions. Plateaus in the boundary
983  indicator (i.e. regions of constant gray value) are handled via a Euclidean distance
984  transform by default.
985 
986  By default, the destination image is assumed to hold seeds for a seeded watershed
987  transform. Seeds may, for example, be created by means of generateWatershedSeeds().
988  Note that the seeds will be overridden with the final watershed segmentation.
989 
990  Alternatively, you may provide \ref SeedOptions in order to instruct
991  watershedsRegionGrowing() to generate its own seeds (it will call generateWatershedSeeds()
992  internally). In that case, the destination image should be zero-initialized.
993 
994  You can specify the neighborhood system to be used by passing \ref FourNeighborCode
995  or \ref EightNeighborCode (default).
996 
997  Further options to be specified via \ref WatershedOptions are:
998 
999  <ul>
1000  <li> Whether to keep a 1-pixel-wide contour (with label 0) between regions or
1001  perform complete grow (i.e. all pixels are assigned to a region).
1002  <li> Whether to stop growing when the boundaryness exceeds a threshold (remaining
1003  pixels keep label 0).
1004  <li> Whether to use a faster, but less powerful algorithm ("turbo algorithm"). It
1005  is faster because it orders pixels by means of a \ref BucketQueue (therefore,
1006  the boundary indicator must contain integers in the range
1007  <tt>[0, ..., bucket_count-1]</tt>, where <tt>bucket_count</tt> is specified in
1008  the options object), it only supports complete growing (no contour between regions
1009  is possible), and it handles plateaus in a simplistic way. It also saves some
1010  memory because it allocates less temporary storage.
1011  <li> Whether one region (label) is to be preferred or discouraged by biasing its cost
1012  with a given factor (smaller than 1 for preference, larger than 1 for discouragement).
1013  </ul>
1014 
1015  Note that VIGRA provides an alternative implementation of the watershed transform via
1016  \ref watershedsUnionFind().
1017 
1018  <b> Declarations:</b>
1019 
1020  pass arguments explicitly:
1021  \code
1022  namespace vigra {
1023  template <class SrcIterator, class SrcAccessor,
1024  class DestIterator, class DestAccessor,
1025  class Neighborhood = EightNeighborCode>
1026  unsigned int
1027  watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1028  DestIterator upperleftd, DestAccessor da,
1029  Neighborhood neighborhood = EightNeighborCode(),
1030  WatershedOptions const & options = WatershedOptions());
1031 
1032  template <class SrcIterator, class SrcAccessor,
1033  class DestIterator, class DestAccessor>
1034  unsigned int
1035  watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1036  DestIterator upperleftd, DestAccessor da,
1037  WatershedOptions const & options = WatershedOptions());
1038  }
1039  \endcode
1040 
1041  use argument objects in conjunction with \ref ArgumentObjectFactories :
1042  \code
1043  namespace vigra {
1044  template <class SrcIterator, class SrcAccessor,
1045  class DestIterator, class DestAccessor,
1046  class Neighborhood = EightNeighborCode>
1047  unsigned int
1048  watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1049  pair<DestIterator, DestAccessor> dest,
1050  Neighborhood neighborhood = EightNeighborCode(),
1051  WatershedOptions const & options = WatershedOptions());
1052 
1053  template <class SrcIterator, class SrcAccessor,
1054  class DestIterator, class DestAccessor>
1055  unsigned int
1056  watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1057  pair<DestIterator, DestAccessor> dest,
1058  WatershedOptions const & options = WatershedOptions());
1059  }
1060  \endcode
1061 
1062  <b> Usage:</b>
1063 
1064  <b>\#include</b> <vigra/watersheds.hxx><br>
1065  Namespace: vigra
1066 
1067  Example: watersheds of the gradient magnitude.
1068 
1069  \code
1070  vigra::BImage src(w, h);
1071  ... // read input data
1072 
1073  // compute gradient magnitude at scale 1.0 as a boundary indicator
1074  vigra::FImage gradMag(w, h);
1075  gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0);
1076 
1077  // example 1
1078  {
1079  // the pixel type of the destination image must be large enough to hold
1080  // numbers up to 'max_region_label' to prevent overflow
1081  vigra::IImage labeling(w, h);
1082 
1083  // call watershed algorithm for 4-neighborhood, leave a 1-pixel boundary between regions,
1084  // and autogenerate seeds from all gradient minima where the magnitude is below 2.0
1085  unsigned int max_region_label =
1086  watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
1087  FourNeighborCode(),
1088  WatershedOptions().keepContours()
1089  .seedOptions(SeedOptions().minima().threshold(2.0)));
1090  }
1091 
1092  // example 2
1093  {
1094  vigra::IImage labeling(w, h);
1095 
1096  // compute seeds beforehand (use connected components of all pixels
1097  // where the gradient is below 4.0)
1098  unsigned int max_region_label =
1099  generateWatershedSeeds(srcImageRange(gradMag), destImage(labeling),
1100  SeedOptions().levelSets(4.0));
1101 
1102  // quantize the gradient image to 256 gray levels
1103  vigra::BImage gradMag256(w, h);
1104  vigra::FindMinMax<float> minmax;
1105  inspectImage(srcImageRange(gradMag), minmax); // find original range
1106  transformImage(srcImageRange(gradMag), destImage(gradMag256),
1107  linearRangeMapping(minmax, 0, 255));
1108 
1109  // call the turbo algorithm with 256 bins, using 8-neighborhood
1110  watershedsRegionGrowing(srcImageRange(gradMag256), destImage(labeling),
1111  WatershedOptions().turboAlgorithm(256));
1112  }
1113 
1114  // example 3
1115  {
1116  vigra::IImage labeling(w, h);
1117 
1118  .. // get seeds from somewhere, e.g. an interactive labeling program,
1119  // make sure that label 1 corresponds to the background
1120 
1121  // bias the watershed algorithm so that the background is preferred
1122  // by reducing the cost for label 1 to 90%
1123  watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
1124  WatershedOptions().biasLabel(1, 0.9));
1125  }
1126  \endcode
1127 
1128  <b> Required Interface:</b>
1129 
1130  \code
1131  SrcIterator src_upperleft, src_lowerright;
1132  DestIterator dest_upperleft;
1133 
1134  SrcAccessor src_accessor;
1135  DestAccessor dest_accessor;
1136 
1137  // compare src values
1138  src_accessor(src_upperleft) <= src_accessor(src_upperleft)
1139 
1140  // set result
1141  int label;
1142  dest_accessor.set(label, dest_upperleft);
1143  \endcode
1144 */
1145 doxygen_overloaded_function(template <...> unsigned int watershedsRegionGrowing)
1146 
1147 template <class SrcIterator, class SrcAccessor,
1148  class DestIterator, class DestAccessor,
1149  class Neighborhood>
1150 unsigned int
1151 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1152  DestIterator upperleftd, DestAccessor da,
1153  Neighborhood neighborhood,
1154  WatershedOptions const & options = WatershedOptions())
1155 {
1156  typedef typename SrcAccessor::value_type ValueType;
1157  typedef typename DestAccessor::value_type LabelType;
1158 
1159  unsigned int max_region_label = 0;
1160 
1161  if(options.seed_options.mini != SeedOptions::Unspecified)
1162  {
1163  // we are supposed to compute seeds
1164  max_region_label =
1165  generateWatershedSeeds(srcIterRange(upperlefts, lowerrights, sa),
1166  destIter(upperleftd, da),
1167  neighborhood, options.seed_options);
1168  }
1169 
1170  if(options.biased_label != 0)
1171  {
1172  // create a statistics functor for biased region growing
1173  detail::BiasedWatershedStatistics<ValueType, LabelType>
1174  regionstats(options.biased_label, options.bias);
1175 
1176  // perform region growing, starting from the seeds computed above
1177  if(options.bucket_count == 0)
1178  {
1179  max_region_label =
1180  seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1181  srcIter(upperleftd, da),
1182  destIter(upperleftd, da),
1183  regionstats, options.terminate, neighborhood, options.max_cost);
1184  }
1185  else
1186  {
1187  max_region_label =
1188  fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1189  destIter(upperleftd, da),
1190  regionstats, options.terminate,
1191  neighborhood, options.max_cost, options.bucket_count);
1192  }
1193  }
1194  else
1195  {
1196  // create a statistics functor for region growing
1197  detail::WatershedStatistics<ValueType, LabelType> regionstats;
1198 
1199  // perform region growing, starting from the seeds computed above
1200  if(options.bucket_count == 0)
1201  {
1202  max_region_label =
1203  seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1204  srcIter(upperleftd, da),
1205  destIter(upperleftd, da),
1206  regionstats, options.terminate, neighborhood, options.max_cost);
1207  }
1208  else
1209  {
1210  max_region_label =
1211  fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1212  destIter(upperleftd, da),
1213  regionstats, options.terminate,
1214  neighborhood, options.max_cost, options.bucket_count);
1215  }
1216  }
1217 
1218  return max_region_label;
1219 }
1220 
1221 template <class SrcIterator, class SrcAccessor,
1222  class DestIterator, class DestAccessor>
1223 inline unsigned int
1224 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1225  DestIterator upperleftd, DestAccessor da,
1226  WatershedOptions const & options = WatershedOptions())
1227 {
1228  return watershedsRegionGrowing(upperlefts, lowerrights, sa, upperleftd, da,
1229  EightNeighborCode(), options);
1230 }
1231 
1232 template <class SrcIterator, class SrcAccessor,
1233  class DestIterator, class DestAccessor,
1234  class Neighborhood>
1235 inline unsigned int
1236 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1237  pair<DestIterator, DestAccessor> dest,
1238  Neighborhood neighborhood,
1239  WatershedOptions const & options = WatershedOptions())
1240 {
1241  return watershedsRegionGrowing(src.first, src.second, src.third,
1242  dest.first, dest.second,
1243  neighborhood, options);
1244 }
1245 
1246 template <class SrcIterator, class SrcAccessor,
1247  class DestIterator, class DestAccessor>
1248 inline unsigned int
1249 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1250  pair<DestIterator, DestAccessor> dest,
1251  WatershedOptions const & options = WatershedOptions())
1252 {
1253  return watershedsRegionGrowing(src.first, src.second, src.third,
1254  dest.first, dest.second,
1255  EightNeighborCode(), options);
1256 }
1257 
1258 
1259 //@}
1260 
1261 } // namespace vigra
1262 
1263 #endif // VIGRA_WATERSHEDS_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (Wed Sep 26 2012)