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

multi_watersheds.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2013 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_MULTI_WATERSHEDS_HXX
37 #define VIGRA_MULTI_WATERSHEDS_HXX
38 
39 #include <functional>
40 #include <limits>
41 #include "mathutil.hxx"
42 #include "multi_array.hxx"
43 #include "multi_math.hxx"
44 #include "multi_gridgraph.hxx"
45 #include "multi_localminmax.hxx"
46 #include "multi_labeling.hxx"
47 #include "watersheds.hxx"
48 #include "bucket_queue.hxx"
49 #include "union_find.hxx"
50 
51 namespace vigra {
52 
53 /** \addtogroup SeededRegionGrowing
54 */
55 //@{
56 namespace lemon_graph {
57 
58 namespace graph_detail {
59 
60  // select the neighbor ID for union-find watersheds
61  // standard behavior: use global node ID
62 template <class Graph>
63 struct NeighborIndexFunctor
64 {
65  typedef typename Graph::index_type index_type;
66 
67  template <class NodeIter, class ArcIter>
68  static index_type get(Graph const & g, NodeIter const &, ArcIter const & a)
69  {
70  return g.id(g.target(*a));
71  }
72 
73  template <class NodeIter, class ArcIter>
74  static index_type getOpposite(Graph const & g, NodeIter const & n, ArcIter const &)
75  {
76  return g.id(*n);
77  }
78 
79  static index_type invalidIndex(Graph const & g)
80  {
81  return std::numeric_limits<index_type>::max();
82  }
83 };
84 
85  // select the neighbor ID for union-find watersheds
86  // GridGraph optimization: use local neighbor index (needs only 1/4 of the memory)
87 template<unsigned int N, class DirectedTag>
88 struct NeighborIndexFunctor<GridGraph<N, DirectedTag> >
89 {
90  typedef GridGraph<N, DirectedTag> Graph;
91  typedef UInt16 index_type;
92 
93  template <class NodeIter, class ArcIter>
94  static index_type get(Graph const & g, NodeIter const &, ArcIter const & a)
95  {
96  return a.neighborIndex();
97  }
98 
99  template <class NodeIter, class ArcIter>
100  static index_type getOpposite(Graph const & g, NodeIter const &, ArcIter const & a)
101  {
102  return g.oppositeIndex(a.neighborIndex());
103  }
104  static index_type invalidIndex(Graph const & g)
105  {
106  return std::numeric_limits<index_type>::max();
107  }
108 };
109 
110 template <class Graph, class T1Map, class T2Map>
111 void
112 prepareWatersheds(Graph const & g,
113  T1Map const & data,
114  T2Map & lowestNeighborIndex)
115 {
116  typedef typename Graph::NodeIt graph_scanner;
117  typedef typename Graph::OutArcIt neighbor_iterator;
118  typedef NeighborIndexFunctor<Graph> IndexFunctor;
119 
120  for (graph_scanner node(g); node != INVALID; ++node)
121  {
122  typename T1Map::value_type lowestValue = data[*node];
123  typename T2Map::value_type lowestIndex = IndexFunctor::invalidIndex(g);
124 
125  for(neighbor_iterator arc(g, *node); arc != INVALID; ++arc)
126  {
127  if(data[g.target(*arc)] < lowestValue)
128  {
129  lowestValue = data[g.target(*arc)];
130  lowestIndex = IndexFunctor::get(g, node, arc);
131  }
132  }
133  lowestNeighborIndex[*node] = lowestIndex;
134  }
135 }
136 
137 
138 template <class Graph, class T1Map, class T2Map, class T3Map>
139 typename T2Map::value_type
140 unionFindWatersheds(Graph const & g,
141  T1Map const & data,
142  T2Map const & lowestNeighborIndex,
143  T3Map & labels)
144 {
145  typedef typename Graph::NodeIt graph_scanner;
146  typedef typename Graph::OutBackArcIt neighbor_iterator;
147  typedef typename T3Map::value_type LabelType;
148  typedef NeighborIndexFunctor<Graph> IndexFunctor;
149 
150  vigra::UnionFindArray<LabelType> regions;
151 
152  // pass 1: find connected components
153  for (graph_scanner node(g); node != INVALID; ++node)
154  {
155  // define tentative label for current node
156  LabelType currentIndex = regions.nextFreeIndex();
157 
158  for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
159  {
160  // merge regions if current target is center's lowest neighbor or vice versa
161  if((lowestNeighborIndex[*node] == IndexFunctor::invalidIndex(g) &&
162  lowestNeighborIndex[g.target(*arc)] == IndexFunctor::invalidIndex(g)) ||
163  (lowestNeighborIndex[*node] == IndexFunctor::get(g, node, arc)) ||
164  (lowestNeighborIndex[g.target(*arc)] == IndexFunctor::getOpposite(g, node, arc)))
165  {
166  currentIndex = regions.makeUnion(labels[g.target(*arc)], currentIndex);
167  }
168  }
169 
170  // set label of current node
171  labels[*node] = regions.finalizeIndex(currentIndex);
172  }
173 
174  LabelType count = regions.makeContiguous();
175 
176  // pass 2: make component labels contiguous
177  for (graph_scanner node(g); node != INVALID; ++node)
178  {
179  labels[*node] = regions.findLabel(labels[*node]);
180  }
181  return count;
182 }
183 
184 template <class Graph, class T1Map, class T2Map>
185 typename T2Map::value_type
186 generateWatershedSeeds(Graph const & g,
187  T1Map const & data,
188  T2Map & seeds,
189  SeedOptions const & options = SeedOptions())
190 {
191  typedef typename T1Map::value_type DataType;
192  typedef unsigned char MarkerType;
193 
194  typename Graph::template NodeMap<MarkerType> minima(g);
195 
196  if(options.mini == SeedOptions::LevelSets)
197  {
198  vigra_precondition(options.thresholdIsValid<DataType>(),
199  "generateWatershedSeeds(): SeedOptions.levelSets() must be specified with threshold.");
200 
201  using namespace multi_math;
202  for(typename Graph::NodeIt iter(g);iter!=lemon::INVALID;++iter){
203  minima[*iter]= data[*iter] <= DataType(options.thresh);
204  }
205  }
206  else
207  {
208  DataType threshold = options.thresholdIsValid<DataType>()
209  ? options.thresh
210  : NumericTraits<DataType>::max();
211 
212  if(options.mini == SeedOptions::ExtendedMinima)
213  extendedLocalMinMaxGraph(g, data, minima, MarkerType(1), threshold,
214  std::less<DataType>(), std::equal_to<DataType>(), true);
215  else
216  lemon_graph::localMinMaxGraph(g, data, minima, MarkerType(1), threshold,
217  std::less<DataType>(), true);
218  }
219  return labelGraphWithBackground(g, minima, seeds, MarkerType(0), std::equal_to<MarkerType>());
220 }
221 
222 
223 template <class Graph, class T1Map, class T2Map>
224 typename T2Map::value_type
225 seededWatersheds(Graph const & g,
226  T1Map const & data,
227  T2Map & labels,
228  WatershedOptions const & options)
229 {
230  typedef typename Graph::Node Node;
231  typedef typename Graph::NodeIt graph_scanner;
232  typedef typename Graph::OutArcIt neighbor_iterator;
233  typedef typename T1Map::value_type CostType;
234  typedef typename T2Map::value_type LabelType;
235 
236  PriorityQueue<Node, CostType, true> pqueue;
237 
238  bool keepContours = ((options.terminate & KeepContours) != 0);
239  LabelType maxRegionLabel = 0;
240 
241  for (graph_scanner node(g); node != INVALID; ++node)
242  {
243  LabelType label = labels[*node];
244  if(label != 0)
245  {
246  if(maxRegionLabel < label)
247  maxRegionLabel = label;
248 
249  for (neighbor_iterator arc(g, *node); arc != INVALID; ++arc)
250  {
251  if(labels[g.target(*arc)] == 0)
252  {
253  // register all seeds that have an unlabeled neighbor
254  if(label == options.biased_label)
255  pqueue.push(*node, data[*node] * options.bias);
256  else
257  pqueue.push(*node, data[*node]);
258  break;
259  }
260  }
261  }
262  }
263 
264  LabelType contourLabel = maxRegionLabel + 1; // temporary contour label
265 
266  // perform region growing
267  while(!pqueue.empty())
268  {
269  Node node = pqueue.top();
270  CostType cost = pqueue.topPriority();
271  pqueue.pop();
272 
273  if((options.terminate & StopAtThreshold) && (cost > options.max_cost))
274  break;
275 
276  LabelType label = labels[node];
277 
278  if(label == contourLabel)
279  continue;
280 
281  // Put the unlabeled neighbors in the priority queue.
282  for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
283  {
284  LabelType neighborLabel = labels[g.target(*arc)];
285  if(neighborLabel == 0)
286  {
287  labels[g.target(*arc)] = label;
288  CostType priority = (label == options.biased_label)
289  ? data[g.target(*arc)] * options.bias
290  : data[g.target(*arc)];
291  if(priority < cost)
292  priority = cost;
293  pqueue.push(g.target(*arc), priority);
294  }
295  else if(keepContours && (label != neighborLabel) && (neighborLabel != contourLabel))
296  {
297  // The present neighbor is adjacent to more than one region
298  // => mark it as contour.
299  CostType priority = (neighborLabel == options.biased_label)
300  ? data[g.target(*arc)] * options.bias
301  : data[g.target(*arc)];
302  if(cost < priority) // neighbor not yet processed
303  labels[g.target(*arc)] = contourLabel;
304  }
305  }
306  }
307 
308  if(keepContours)
309  {
310  // Replace the temporary contour label with label 0.
311  ///typename T2Map::iterator k = labels.begin(),
312  /// end = labels.end();
313  ///for(; k != end; ++k)
314  /// if(*k == contourLabel)
315  /// *k = 0;
316 
317  for(typename Graph::NodeIt iter(g);iter!=lemon::INVALID;++iter){
318  if(labels[*iter]==contourLabel)
319  labels[*iter]=0;
320  }
321  }
322 
323  return maxRegionLabel;
324 }
325 
326 } // namespace graph_detail
327 
328 template <class Graph, class T1Map, class T2Map>
329 typename T2Map::value_type
330 watershedsGraph(Graph const & g,
331  T1Map const & data,
332  T2Map & labels,
333  WatershedOptions const & options)
334 {
335  if(options.method == WatershedOptions::UnionFind)
336  {
337  typedef typename graph_detail::NeighborIndexFunctor<Graph>::index_type index_type;
338 
339  vigra_precondition((index_type)g.maxDegree() <= NumericTraits<index_type>::max(),
340  "watershedsGraph(): cannot handle nodes with degree > 65535.");
341 
342  typename Graph::template NodeMap<index_type> lowestNeighborIndex(g);
343 
344  graph_detail::prepareWatersheds(g, data, lowestNeighborIndex);
345  return graph_detail::unionFindWatersheds(g, data, lowestNeighborIndex, labels);
346  }
347  else if(options.method == WatershedOptions::RegionGrowing)
348  {
349  SeedOptions seed_options;
350 
351  // check if the user has explicitly requested seed computation
352  if(options.seed_options.mini != SeedOptions::Unspecified)
353  {
354  seed_options = options.seed_options;
355  }
356  else
357  {
358  // otherwise, don't compute seeds if 'labels' already contains them
359  if(labels.any())
360  seed_options.mini = SeedOptions::Unspecified;
361  }
362 
363  if(seed_options.mini != SeedOptions::Unspecified)
364  {
365  graph_detail::generateWatershedSeeds(g, data, labels, seed_options);
366  }
367 
368  return graph_detail::seededWatersheds(g, data, labels, options);
369  }
370  else
371  {
372  vigra_precondition(false,
373  "watershedsGraph(): invalid method in watershed options.");
374  return 0;
375  }
376 }
377 
378 
379 } // namespace lemon_graph
380 
381  // documentation is in watersheds.hxx
382 template <unsigned int N, class T, class S1,
383  class Label, class S2>
384 inline Label
385 generateWatershedSeeds(MultiArrayView<N, T, S1> const & data,
386  MultiArrayView<N, Label, S2> seeds,
387  NeighborhoodType neighborhood = IndirectNeighborhood,
388  SeedOptions const & options = SeedOptions())
389 {
390  vigra_precondition(data.shape() == seeds.shape(),
391  "generateWatershedSeeds(): Shape mismatch between input and output.");
392 
393  GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
394  return lemon_graph::graph_detail::generateWatershedSeeds(graph, data, seeds, options);
395 }
396 
397 
398 /** \brief Watershed segmentation of an arbitrary-dimensional array.
399 
400  See also \ref unionFindWatershedsBlockwise() for a parallel version of the
401  watershed algorithm.
402 
403  This function implements variants of the watershed algorithms
404  described in
405 
406  [1] L. Vincent and P. Soille: <em>"Watersheds in digital spaces: An efficient algorithm
407  based on immersion simulations"</em>, IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991
408 
409  [2] J. Roerdink, R. Meijster: <em>"The watershed transform: definitions, algorithms,
410  and parallelization strategies"</em>, Fundamenta Informaticae, 41:187-228, 2000
411 
412  The source array \a data is a boundary indicator such as the gaussianGradientMagnitude()
413  or the trace of the \ref boundaryTensor(), and the destination \a labels is a label array
414  designating membership of each point in one of the regions found. Plateaus in the boundary
415  indicator are handled via simple tie breaking strategies. Argument \a neighborhood
416  specifies the connectivity between points and can be <tt>DirectNeighborhood</tt> (meaning
417  4-neighborhood in 2D and 6-neighborhood in 3D, default) or <tt>IndirectNeighborhood</tt>
418  (meaning 8-neighborhood in 2D and 26-neighborhood in 3D).
419 
420  The watershed variant to be applied can be selected in the \ref WatershedOptions
421  object: When you call <tt>WatershedOptions::regionGrowing()</tt> (default), the flooding
422  algorithm from [1] is used. Alternatively, <tt>WatershedOptions::unionFind()</tt> uses
423  the scan-line algorithm 4.7 from [2]. The latter is faster, but does not support any options
424  (if you pass options nonetheless, they are silently ignored).
425 
426  The region growing algorithm needs a seed for each region. Seeds can either be provided in
427  the destination array \a labels (which will then be overwritten with the result) or computed
428  automatically by an internal call to generateWatershedSeeds(). In the former case you have
429  full control over seed placement, while the latter is more convenient. Automatic seed
430  computation is performed when you provide seeding options via <tt>WatershedOptions::seedOptions()</tt>
431  or when the array \a labels is empty (all zeros), in which case default seeding options
432  are chosen. The destination image should be zero-initialized for automatic seed computation.
433 
434  Further options to be specified via \ref WatershedOptions are:
435 
436  <ul>
437  <li> <tt>keepContours()</tt>: Whether to keep a 1-pixel-wide contour (with label 0) between
438  regions (otherwise, a complete grow is performed, i.e. all pixels are assigned to a region).
439  <li> <tt>stopAtThreshold()</tt>: Whether to stop growing when the boundaryness exceeds a threshold
440  (remaining pixels keep label 0).
441  <li> <tt>biasLabel()</tt>: Whether one region (label) is to be preferred or discouraged by biasing its cost
442  with a given factor (smaller than 1 for preference, larger than 1 for discouragement).
443  </ul>
444 
445  The option <tt>turboAlgorithm()</tt> is implied by method <tt>regionGrowing()</tt> (this is
446  in contrast to watershedsRegionGrowing(), which supports an additional algorithm in 2D only).
447 
448  watershedsMultiArray() returns the number of regions found (= the highest region label, because
449  labels start at 1).
450 
451  <b> Declaration:</b>
452 
453  \code
454  namespace vigra {
455  template <unsigned int N, class T, class S1,
456  class Label, class S2>
457  Label
458  watershedsMultiArray(MultiArrayView<N, T, S1> const & data,
459  MultiArrayView<N, Label, S2> labels, // may also hold input seeds
460  NeighborhoodType neighborhood = DirectNeighborhood,
461  WatershedOptions const & options = WatershedOptions());
462  }
463  \endcode
464 
465  <b> Usage:</b>
466 
467  <b>\#include</b> <vigra/multi_watersheds.hxx><br>
468  Namespace: vigra
469 
470  Example: watersheds of the gradient magnitude (the example works likewise for higher dimensions).
471 
472  \code
473  MultiArray<2, unsigned char> src(Shape2(w, h));
474  ... // read input data
475 
476  // compute gradient magnitude at scale 1.0 as a boundary indicator
477  MultiArray<2, float> gradMag(src.shape());
478  gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0);
479 
480  // example 1
481  {
482  // the pixel type of the destination image must be large enough to hold
483  // numbers up to 'max_region_label' to prevent overflow
484  MultiArray<2, unsigned int> labeling(src.shape());
485 
486  // call region-growing algorithm for 4-neighborhood, leave a 1-pixel boundary between
487  // regions, and autogenerate seeds from all gradient minima where the magnitude is
488  // less than 2.0.
489  unsigned int max_region_label =
490  watershedsMultiArray(gradMag, labeling, DirectNeighborhood,
491  WatershedOptions().keepContours()
492  .seedOptions(SeedOptions().minima().threshold(2.0)));
493  }
494 
495  // example 2
496  {
497  MultiArray<2, unsigned int> labeling(src.shape());
498 
499  // compute seeds beforehand (use connected components of all pixels
500  // where the gradient is below 4.0)
501  unsigned int max_region_label = generateWatershedSeeds(gradMag, labeling,
502  SeedOptions().levelSets(4.0));
503 
504  // quantize the gradient image to 256 gray levels
505  float m, M;
506  gradMag.minmax(&m, &M);
507 
508  using namespace multi_math;
509  MultiArray<2, unsigned char> gradMag256(255.0 / (M - m) * (gradMag - m));
510 
511  // call region-growing algorithm with 8-neighborhood,
512  // since the data are 8-bit, a faster priority queue will be used
513  watershedsMultiArray(gradMag256, labeling, IndirectNeighborhood);
514  }
515 
516  // example 3
517  {
518  MultiArray<2, unsigned int> labeling(src.shape());
519 
520  .. // put seeds in 'labeling', e.g. from an interactive labeling program,
521  // make sure that label 1 corresponds to the background
522 
523  // bias the watershed algorithm so that the background is preferred
524  // by reducing the cost for label 1 to 90%
525  watershedsMultiArray(gradMag, labeling,
526  WatershedOptions().biasLabel(1, 0.9));
527  }
528 
529  // example 4
530  {
531  MultiArray<2, unsigned int> labeling(src.shape());
532 
533  // use the fast union-find algorithm with 4-neighborhood
534  watershedsMultiArray(gradMag, labeling, WatershedOptions().unionFind());
535  }
536  \endcode
537 */
539 
540 template <unsigned int N, class T, class S1,
541  class Label, class S2>
542 inline Label
543 watershedsMultiArray(MultiArrayView<N, T, S1> const & data,
544  MultiArrayView<N, Label, S2> labels, // may also hold input seeds
545  NeighborhoodType neighborhood = DirectNeighborhood,
546  WatershedOptions const & options = WatershedOptions())
547 {
548  vigra_precondition(data.shape() == labels.shape(),
549  "watershedsMultiArray(): Shape mismatch between input and output.");
550 
551  GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
552  return lemon_graph::watershedsGraph(graph, data, labels, options);
553 }
554 
555 //@}
556 
557 } // namespace vigra
558 
559 #endif // VIGRA_MULTI_WATERSHEDS_HXX
CoupledHandleCast< TARGET_INDEX, Handle >::reference get(Handle &handle)
Definition: multi_handle.hxx:927
detail::SelectIntegerType< 16, detail::UnsignedIntTypes >::type UInt16
16-bit unsigned int
Definition: sized_int.hxx:181
Definition: accessor.hxx:43
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
Label watershedsMultiArray(...)
Watershed segmentation of an arbitrary-dimensional array.
use direct and indirect neighbors
Definition: multi_fwd.hxx:188
use only direct neighbors
Definition: multi_fwd.hxx:187
NeighborhoodType
Choose the neighborhood system in a dimension-independent way.
Definition: multi_fwd.hxx:186

© 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.10.0