[ 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  {
405  return thresh < double(NumericTraits<T>::max());
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  LocalMinmaxOptions lm_options;
505  lm_options.neighborhood(Neighborhood::DirectionCount)
506  .markWith(1.0)
507  .allowAtBorder()
508  .allowPlateaus(options.mini == SeedOptions::ExtendedMinima);
509  if(options.thresholdIsValid<SrcType>())
510  lm_options.threshold(options.thresh);
511 
512  localMinima(srcIterRange(upperlefts, lowerrights, sa), destImage(seeds),
513  lm_options);
514  }
515 
516  return labelImageWithBackground(srcImageRange(seeds), destIter(upperleftd, da),
517  Neighborhood::DirectionCount == 8, 0);
518 }
519 
520 template <class SrcIterator, class SrcAccessor,
521  class DestIterator, class DestAccessor>
522 inline unsigned int
523 generateWatershedSeeds(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
524  DestIterator upperleftd, DestAccessor da,
525  SeedOptions const & options = SeedOptions())
526 {
527  return generateWatershedSeeds(upperlefts, lowerrights, sa, upperleftd, da,
528  EightNeighborCode(), options);
529 }
530 
531 template <class SrcIterator, class SrcAccessor,
532  class DestIterator, class DestAccessor,
533  class Neighborhood>
534 inline unsigned int
535 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
536  pair<DestIterator, DestAccessor> dest,
537  Neighborhood neighborhood,
538  SeedOptions const & options = SeedOptions())
539 {
540  return generateWatershedSeeds(src.first, src.second, src.third,
541  dest.first, dest.second,
542  neighborhood, options);
543 }
544 
545 template <class SrcIterator, class SrcAccessor,
546  class DestIterator, class DestAccessor>
547 inline unsigned int
548 generateWatershedSeeds(triple<SrcIterator, SrcIterator, SrcAccessor> src,
549  pair<DestIterator, DestAccessor> dest,
550  SeedOptions const & options = SeedOptions())
551 {
552  return generateWatershedSeeds(src.first, src.second, src.third,
553  dest.first, dest.second,
554  EightNeighborCode(), options);
555 }
556 
557 /********************************************************/
558 /* */
559 /* watersheds */
560 /* */
561 /********************************************************/
562 
563 /** \brief Region segmentation by means of the union-find watershed algorithm.
564 
565  This function implements the union-find version of the watershed algorithms
566  as described in
567 
568  J. Roerdink, R. Meijster: "<em>The watershed transform: definitions, algorithms,
569  and parallelization strategies</em>", Fundamenta Informaticae, 41:187-228, 2000
570 
571  The source image is a boundary indicator such as the gaussianGradientMagnitude()
572  or the trace of the \ref boundaryTensor(). Local minima of the boundary indicator
573  are used as region seeds, and all other pixels are recursively assigned to the same
574  region as their lowest neighbor. Pass \ref vigra::EightNeighborCode or
575  \ref vigra::FourNeighborCode to determine the neighborhood where pixel values
576  are compared. The pixel type of the input image must be <tt>LessThanComparable</tt>.
577  The function uses accessors.
578 
579  Note that VIGRA provides an alternative implementation of the watershed transform via
580  \ref watershedsRegionGrowing(). It is slower, but offers many more configuration options.
581 
582  <b> Declarations:</b>
583 
584  pass arguments explicitly:
585  \code
586  namespace vigra {
587  template <class SrcIterator, class SrcAccessor,
588  class DestIterator, class DestAccessor,
589  class Neighborhood = EightNeighborCode>
590  unsigned int
591  watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
592  DestIterator upperleftd, DestAccessor da,
593  Neighborhood neighborhood = EightNeighborCode())
594  }
595  \endcode
596 
597  use argument objects in conjunction with \ref ArgumentObjectFactories :
598  \code
599  namespace vigra {
600  template <class SrcIterator, class SrcAccessor,
601  class DestIterator, class DestAccessor,
602  class Neighborhood = EightNeighborCode>
603  unsigned int
604  watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
605  pair<DestIterator, DestAccessor> dest,
606  Neighborhood neighborhood = EightNeighborCode())
607  }
608  \endcode
609 
610  <b> Usage:</b>
611 
612  <b>\#include</b> <vigra/watersheds.hxx><br>
613  Namespace: vigra
614 
615  Example: watersheds of the gradient magnitude.
616 
617  \code
618  vigra::BImage in(w,h);
619  ... // read input data
620 
621  // compute gradient magnitude as boundary indicator
622  vigra::FImage gradMag(w, h);
623  gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 3.0);
624 
625  // the pixel type of the destination image must be large enough to hold
626  // numbers up to 'max_region_label' to prevent overflow
627  vigra::IImage labeling(w,h);
628  int max_region_label = watershedsUnionFind(srcImageRange(gradMag), destImage(labeling));
629 
630  \endcode
631 
632  <b> Required Interface:</b>
633 
634  \code
635  SrcIterator src_upperleft, src_lowerright;
636  DestIterator dest_upperleft;
637 
638  SrcAccessor src_accessor;
639  DestAccessor dest_accessor;
640 
641  // compare src values
642  src_accessor(src_upperleft) <= src_accessor(src_upperleft)
643 
644  // set result
645  int label;
646  dest_accessor.set(label, dest_upperleft);
647  \endcode
648 */
649 doxygen_overloaded_function(template <...> unsigned int watershedsUnionFind)
650 
651 template <class SrcIterator, class SrcAccessor,
652  class DestIterator, class DestAccessor,
653  class Neighborhood>
654 unsigned int
655 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
656  DestIterator upperleftd, DestAccessor da,
657  Neighborhood neighborhood)
658 {
659  SImage orientationImage(lowerrights - upperlefts);
660 
661  prepareWatersheds(upperlefts, lowerrights, sa,
662  orientationImage.upperLeft(), orientationImage.accessor(), neighborhood);
663  return watershedLabeling(orientationImage.upperLeft(), orientationImage.lowerRight(), orientationImage.accessor(),
664  upperleftd, da, neighborhood);
665 }
666 
667 template <class SrcIterator, class SrcAccessor,
668  class DestIterator, class DestAccessor>
669 inline unsigned int
670 watershedsUnionFind(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
671  DestIterator upperleftd, DestAccessor da)
672 {
673  return watershedsUnionFind(upperlefts, lowerrights, sa, upperleftd, da, EightNeighborCode());
674 }
675 
676 template <class SrcIterator, class SrcAccessor,
677  class DestIterator, class DestAccessor,
678  class Neighborhood>
679 inline unsigned int
680 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
681  pair<DestIterator, DestAccessor> dest, Neighborhood neighborhood)
682 {
683  return watershedsUnionFind(src.first, src.second, src.third,
684  dest.first, dest.second, neighborhood);
685 }
686 
687 template <class SrcIterator, class SrcAccessor,
688  class DestIterator, class DestAccessor>
689 inline unsigned int
690 watershedsUnionFind(triple<SrcIterator, SrcIterator, SrcAccessor> src,
691  pair<DestIterator, DestAccessor> dest)
692 {
693  return watershedsUnionFind(src.first, src.second, src.third,
694  dest.first, dest.second);
695 }
696 
697 /** \brief Options object for watershedsRegionGrowing().
698 
699  <b> Usage:</b>
700 
701  see watershedsRegionGrowing() for detailed examples.
702 */
704 {
705  public:
706  double max_cost, bias;
707  SRGType terminate;
708  unsigned int biased_label, bucket_count;
709  SeedOptions seed_options;
710 
711  /** \brief Create options object with default settings.
712 
713  Defaults are: perform complete grow (all pixels are assigned to regions),
714  use standard algorithm, assume that the destination image already contains
715  region seeds.
716  */
718  : max_cost(0.0),
719  bias(1.0),
720  terminate(CompleteGrow),
721  biased_label(0),
722  bucket_count(0),
723  seed_options(SeedOptions().unspecified())
724  {}
725 
726  /** \brief Perform complete grow.
727 
728  That is, all pixels are assigned to regions, without explicit contours
729  in between.
730 
731  Default: true
732  */
734  {
735  terminate = SRGType(CompleteGrow | (terminate & StopAtThreshold));
736  return *this;
737  }
738 
739  /** \brief Keep one-pixel wide contour between regions.
740 
741  Note that this option is unsupported by the turbo algorithm.
742 
743  Default: false
744  */
746  {
747  terminate = SRGType(KeepContours | (terminate & StopAtThreshold));
748  return *this;
749  }
750 
751  /** \brief Set \ref SRGType explicitly.
752 
753  Default: CompleteGrow
754  */
756  {
757  terminate = type;
758  return *this;
759  }
760 
761  /** \brief Stop region growing when the boundaryness exceeds the threshold.
762 
763  This option may be combined with completeGrow() and keepContours().
764 
765  Default: no early stopping
766  */
767  WatershedOptions & stopAtThreshold(double threshold)
768  {
769  terminate = SRGType(terminate | StopAtThreshold);
770  max_cost = threshold;
771  return *this;
772  }
773 
774  /** \brief Use a simpler, but faster region growing algorithm.
775 
776  The algorithm internally uses a \ref BucketQueue to determine
777  the processing order of the pixels. This is only useful,
778  when the input boundary indicator image contains integers
779  in the range <tt>[0, ..., bucket_count-1]</tt>. Since
780  these boundary indicators are typically represented as
781  UInt8 images, the default <tt>bucket_count</tt> is 256.
782 
783  Default: don't use the turbo algorithm
784  */
785  WatershedOptions & turboAlgorithm(unsigned int bucket_count = 256)
786  {
787  this->bucket_count = bucket_count;
788  return *this;
789  }
790 
791  /** \brief Specify seed options.
792 
793  In this case, watershedsRegionGrowing() assumes that the destination
794  image does not yet contain seeds. It will therefore call
795  generateWatershedSeeds() and pass on the seed options.
796 
797  Default: don't compute seeds (i.e. assume that destination image already
798  contains seeds).
799  */
801  {
802  seed_options = s;
803  return *this;
804  }
805 
806  /** \brief Bias the cost of the specified region by the given factor.
807 
808  In certain applications, one region (typically the background) should
809  be preferred in region growing. This is most easily achieved
810  by adjusting the assignment cost for that region as <tt>factor*cost</tt>,
811  with a factor slightly below 1.
812 
813  Default: don't bias any region.
814  */
815  WatershedOptions & biasLabel(unsigned int label, double factor)
816  {
817  biased_label = label;
818  bias = factor;
819  return *this;
820  }
821 };
822 
823 namespace detail {
824 
825 template <class CostType, class LabelType>
826 class WatershedStatistics
827 {
828  public:
829 
830  typedef SeedRgDirectValueFunctor<CostType> value_type;
831  typedef value_type & reference;
832  typedef value_type const & const_reference;
833 
834  typedef CostType first_argument_type;
835  typedef LabelType second_argument_type;
836  typedef LabelType argument_type;
837 
838  WatershedStatistics()
839  {}
840 
841  void resize(unsigned int)
842  {}
843 
844  void reset()
845  {}
846 
847  /** update regions statistics (do nothing in the watershed algorithm)
848  */
849  template <class T1, class T2>
850  void operator()(first_argument_type const &, second_argument_type const &)
851  {}
852 
853  /** ask for maximal index (label) allowed
854  */
855  LabelType maxRegionLabel() const
856  { return size() - 1; }
857 
858  /** ask for array size (i.e. maxRegionLabel() + 1)
859  */
860  LabelType size() const
861  { return NumericTraits<LabelType>::max(); }
862 
863  /** read the statistics functor for a region via its label
864  */
865  const_reference operator[](argument_type label) const
866  { return stats; }
867 
868  /** access the statistics functor for a region via its label
869  */
870  reference operator[](argument_type label)
871  { return stats; }
872 
873  value_type stats;
874 };
875 
876 template <class Value>
877 class SeedRgBiasedValueFunctor
878 {
879  public:
880  double bias;
881 
882  /* the functor's argument type
883  */
884  typedef Value argument_type;
885 
886  /* the functor's result type (unused, only necessary for
887  use of SeedRgDirectValueFunctor in \ref vigra::ArrayOfRegionStatistics
888  */
889  typedef Value result_type;
890 
891  /* the return type of the cost() function
892  */
893  typedef Value cost_type;
894 
895  SeedRgBiasedValueFunctor(double b = 1.0)
896  : bias(b)
897  {}
898 
899  /* Do nothing (since we need not update region statistics).
900  */
901  void operator()(argument_type const &) const {}
902 
903  /* Return scaled argument
904  */
905  cost_type cost(argument_type const & v) const
906  {
907  return cost_type(bias*v);
908  }
909 };
910 
911 template <class CostType, class LabelType>
912 class BiasedWatershedStatistics
913 {
914  public:
915 
916  typedef SeedRgBiasedValueFunctor<CostType> value_type;
917  typedef value_type & reference;
918  typedef value_type const & const_reference;
919 
920  typedef CostType first_argument_type;
921  typedef LabelType second_argument_type;
922  typedef LabelType argument_type;
923 
924  BiasedWatershedStatistics(LabelType biasedLabel, double bias)
925  : biased_label(biasedLabel),
926  biased_stats(bias)
927  {}
928 
929  void resize(unsigned int)
930  {}
931 
932  void reset()
933  {}
934 
935  /** update regions statistics (do nothing in the watershed algorithm)
936  */
937  template <class T1, class T2>
938  void operator()(first_argument_type const &, second_argument_type const &)
939  {}
940 
941  /** ask for maximal index (label) allowed
942  */
943  LabelType maxRegionLabel() const
944  { return size() - 1; }
945 
946  /** ask for array size (i.e. maxRegionLabel() + 1)
947  */
948  LabelType size() const
949  { return NumericTraits<LabelType>::max(); }
950 
951  /** read the statistics functor for a region via its label
952  */
953  const_reference operator[](argument_type label) const
954  {
955  return (label == biased_label)
956  ? biased_stats
957  : stats;
958  }
959 
960  /** access the statistics functor for a region via its label
961  */
962  reference operator[](argument_type label)
963  {
964  return (label == biased_label)
965  ? biased_stats
966  : stats;
967  }
968 
969  LabelType biased_label;
970  value_type stats, biased_stats;
971 };
972 
973 } // namespace detail
974 
975 /** \brief Region segmentation by means of a flooding-based watershed algorithm.
976 
977  This function implements variants of the watershed algorithm
978  described in
979 
980  L. Vincent and P. Soille: "<em>Watersheds in digital spaces: An efficient algorithm
981  based on immersion simulations</em>", IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991
982 
983  The source image is a boundary indicator such as the gaussianGradientMagnitude()
984  or the trace of the \ref boundaryTensor(), and the destination is a label image
985  designating membership of each pixel in one of the regions. Plateaus in the boundary
986  indicator (i.e. regions of constant gray value) are handled via a Euclidean distance
987  transform by default.
988 
989  By default, the destination image is assumed to hold seeds for a seeded watershed
990  transform. Seeds may, for example, be created by means of generateWatershedSeeds().
991  Note that the seeds will be overridden with the final watershed segmentation.
992 
993  Alternatively, you may provide \ref SeedOptions in order to instruct
994  watershedsRegionGrowing() to generate its own seeds (it will call generateWatershedSeeds()
995  internally). In that case, the destination image should be zero-initialized.
996 
997  You can specify the neighborhood system to be used by passing \ref FourNeighborCode
998  or \ref EightNeighborCode (default).
999 
1000  Further options to be specified via \ref WatershedOptions are:
1001 
1002  <ul>
1003  <li> Whether to keep a 1-pixel-wide contour (with label 0) between regions or
1004  perform complete grow (i.e. all pixels are assigned to a region).
1005  <li> Whether to stop growing when the boundaryness exceeds a threshold (remaining
1006  pixels keep label 0).
1007  <li> Whether to use a faster, but less powerful algorithm ("turbo algorithm"). It
1008  is faster because it orders pixels by means of a \ref BucketQueue (therefore,
1009  the boundary indicator must contain integers in the range
1010  <tt>[0, ..., bucket_count-1]</tt>, where <tt>bucket_count</tt> is specified in
1011  the options object), it only supports complete growing (no contour between regions
1012  is possible), and it handles plateaus in a simplistic way. It also saves some
1013  memory because it allocates less temporary storage.
1014  <li> Whether one region (label) is to be preferred or discouraged by biasing its cost
1015  with a given factor (smaller than 1 for preference, larger than 1 for discouragement).
1016  </ul>
1017 
1018  Note that VIGRA provides an alternative implementation of the watershed transform via
1019  \ref watershedsUnionFind().
1020 
1021  <b> Declarations:</b>
1022 
1023  pass arguments explicitly:
1024  \code
1025  namespace vigra {
1026  template <class SrcIterator, class SrcAccessor,
1027  class DestIterator, class DestAccessor,
1028  class Neighborhood = EightNeighborCode>
1029  unsigned int
1030  watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1031  DestIterator upperleftd, DestAccessor da,
1032  Neighborhood neighborhood = EightNeighborCode(),
1033  WatershedOptions const & options = WatershedOptions());
1034 
1035  template <class SrcIterator, class SrcAccessor,
1036  class DestIterator, class DestAccessor>
1037  unsigned int
1038  watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1039  DestIterator upperleftd, DestAccessor da,
1040  WatershedOptions const & options = WatershedOptions());
1041  }
1042  \endcode
1043 
1044  use argument objects in conjunction with \ref ArgumentObjectFactories :
1045  \code
1046  namespace vigra {
1047  template <class SrcIterator, class SrcAccessor,
1048  class DestIterator, class DestAccessor,
1049  class Neighborhood = EightNeighborCode>
1050  unsigned int
1051  watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1052  pair<DestIterator, DestAccessor> dest,
1053  Neighborhood neighborhood = EightNeighborCode(),
1054  WatershedOptions const & options = WatershedOptions());
1055 
1056  template <class SrcIterator, class SrcAccessor,
1057  class DestIterator, class DestAccessor>
1058  unsigned int
1059  watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1060  pair<DestIterator, DestAccessor> dest,
1061  WatershedOptions const & options = WatershedOptions());
1062  }
1063  \endcode
1064 
1065  <b> Usage:</b>
1066 
1067  <b>\#include</b> <vigra/watersheds.hxx><br>
1068  Namespace: vigra
1069 
1070  Example: watersheds of the gradient magnitude.
1071 
1072  \code
1073  vigra::BImage src(w, h);
1074  ... // read input data
1075 
1076  // compute gradient magnitude at scale 1.0 as a boundary indicator
1077  vigra::FImage gradMag(w, h);
1078  gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0);
1079 
1080  // example 1
1081  {
1082  // the pixel type of the destination image must be large enough to hold
1083  // numbers up to 'max_region_label' to prevent overflow
1084  vigra::IImage labeling(w, h);
1085 
1086  // call watershed algorithm for 4-neighborhood, leave a 1-pixel boundary between regions,
1087  // and autogenerate seeds from all gradient minima where the magnitude is below 2.0
1088  unsigned int max_region_label =
1089  watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
1090  FourNeighborCode(),
1091  WatershedOptions().keepContours()
1092  .seedOptions(SeedOptions().minima().threshold(2.0)));
1093  }
1094 
1095  // example 2
1096  {
1097  vigra::IImage labeling(w, h);
1098 
1099  // compute seeds beforehand (use connected components of all pixels
1100  // where the gradient is below 4.0)
1101  unsigned int max_region_label =
1102  generateWatershedSeeds(srcImageRange(gradMag), destImage(labeling),
1103  SeedOptions().levelSets(4.0));
1104 
1105  // quantize the gradient image to 256 gray levels
1106  vigra::BImage gradMag256(w, h);
1107  vigra::FindMinMax<float> minmax;
1108  inspectImage(srcImageRange(gradMag), minmax); // find original range
1109  transformImage(srcImageRange(gradMag), destImage(gradMag256),
1110  linearRangeMapping(minmax, 0, 255));
1111 
1112  // call the turbo algorithm with 256 bins, using 8-neighborhood
1113  watershedsRegionGrowing(srcImageRange(gradMag256), destImage(labeling),
1114  WatershedOptions().turboAlgorithm(256));
1115  }
1116 
1117  // example 3
1118  {
1119  vigra::IImage labeling(w, h);
1120 
1121  .. // get seeds from somewhere, e.g. an interactive labeling program,
1122  // make sure that label 1 corresponds to the background
1123 
1124  // bias the watershed algorithm so that the background is preferred
1125  // by reducing the cost for label 1 to 90%
1126  watershedsRegionGrowing(srcImageRange(gradMag), destImage(labeling),
1127  WatershedOptions().biasLabel(1, 0.9));
1128  }
1129  \endcode
1130 
1131  <b> Required Interface:</b>
1132 
1133  \code
1134  SrcIterator src_upperleft, src_lowerright;
1135  DestIterator dest_upperleft;
1136 
1137  SrcAccessor src_accessor;
1138  DestAccessor dest_accessor;
1139 
1140  // compare src values
1141  src_accessor(src_upperleft) <= src_accessor(src_upperleft)
1142 
1143  // set result
1144  int label;
1145  dest_accessor.set(label, dest_upperleft);
1146  \endcode
1147 */
1148 doxygen_overloaded_function(template <...> unsigned int watershedsRegionGrowing)
1149 
1150 template <class SrcIterator, class SrcAccessor,
1151  class DestIterator, class DestAccessor,
1152  class Neighborhood>
1153 unsigned int
1154 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1155  DestIterator upperleftd, DestAccessor da,
1156  Neighborhood neighborhood,
1157  WatershedOptions const & options = WatershedOptions())
1158 {
1159  typedef typename SrcAccessor::value_type ValueType;
1160  typedef typename DestAccessor::value_type LabelType;
1161 
1162  unsigned int max_region_label = 0;
1163 
1164  if(options.seed_options.mini != SeedOptions::Unspecified)
1165  {
1166  // we are supposed to compute seeds
1167  max_region_label =
1168  generateWatershedSeeds(srcIterRange(upperlefts, lowerrights, sa),
1169  destIter(upperleftd, da),
1170  neighborhood, options.seed_options);
1171  }
1172 
1173  if(options.biased_label != 0)
1174  {
1175  // create a statistics functor for biased region growing
1176  detail::BiasedWatershedStatistics<ValueType, LabelType>
1177  regionstats(options.biased_label, options.bias);
1178 
1179  // perform region growing, starting from the seeds computed above
1180  if(options.bucket_count == 0)
1181  {
1182  max_region_label =
1183  seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1184  srcIter(upperleftd, da),
1185  destIter(upperleftd, da),
1186  regionstats, options.terminate, neighborhood, options.max_cost);
1187  }
1188  else
1189  {
1190  max_region_label =
1191  fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1192  destIter(upperleftd, da),
1193  regionstats, options.terminate,
1194  neighborhood, options.max_cost, options.bucket_count);
1195  }
1196  }
1197  else
1198  {
1199  // create a statistics functor for region growing
1200  detail::WatershedStatistics<ValueType, LabelType> regionstats;
1201 
1202  // perform region growing, starting from the seeds computed above
1203  if(options.bucket_count == 0)
1204  {
1205  max_region_label =
1206  seededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1207  srcIter(upperleftd, da),
1208  destIter(upperleftd, da),
1209  regionstats, options.terminate, neighborhood, options.max_cost);
1210  }
1211  else
1212  {
1213  max_region_label =
1214  fastSeededRegionGrowing(srcIterRange(upperlefts, lowerrights, sa),
1215  destIter(upperleftd, da),
1216  regionstats, options.terminate,
1217  neighborhood, options.max_cost, options.bucket_count);
1218  }
1219  }
1220 
1221  return max_region_label;
1222 }
1223 
1224 template <class SrcIterator, class SrcAccessor,
1225  class DestIterator, class DestAccessor>
1226 inline unsigned int
1227 watershedsRegionGrowing(SrcIterator upperlefts, SrcIterator lowerrights, SrcAccessor sa,
1228  DestIterator upperleftd, DestAccessor da,
1229  WatershedOptions const & options = WatershedOptions())
1230 {
1231  return watershedsRegionGrowing(upperlefts, lowerrights, sa, upperleftd, da,
1232  EightNeighborCode(), options);
1233 }
1234 
1235 template <class SrcIterator, class SrcAccessor,
1236  class DestIterator, class DestAccessor,
1237  class Neighborhood>
1238 inline unsigned int
1239 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1240  pair<DestIterator, DestAccessor> dest,
1241  Neighborhood neighborhood,
1242  WatershedOptions const & options = WatershedOptions())
1243 {
1244  return watershedsRegionGrowing(src.first, src.second, src.third,
1245  dest.first, dest.second,
1246  neighborhood, options);
1247 }
1248 
1249 template <class SrcIterator, class SrcAccessor,
1250  class DestIterator, class DestAccessor>
1251 inline unsigned int
1252 watershedsRegionGrowing(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1253  pair<DestIterator, DestAccessor> dest,
1254  WatershedOptions const & options = WatershedOptions())
1255 {
1256  return watershedsRegionGrowing(src.first, src.second, src.third,
1257  dest.first, dest.second,
1258  EightNeighborCode(), options);
1259 }
1260 
1261 
1262 //@}
1263 
1264 } // namespace vigra
1265 
1266 #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.9.0 (Wed Feb 27 2013)