[ 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  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(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  This function implements variants of the watershed algorithms
401  described in
402 
403  [1] L. Vincent and P. Soille: <em>"Watersheds in digital spaces: An efficient algorithm
404  based on immersion simulations"</em>, IEEE Trans. Patt. Analysis Mach. Intell. 13(6):583-598, 1991
405 
406  [2] J. Roerdink, R. Meijster: <em>"The watershed transform: definitions, algorithms,
407  and parallelization strategies"</em>, Fundamenta Informaticae, 41:187-228, 2000
408 
409  The source array \a data is a boundary indicator such as the gaussianGradientMagnitude()
410  or the trace of the \ref boundaryTensor(), and the destination \a labels is a label array
411  designating membership of each point in one of the regions found. Plateaus in the boundary
412  indicator are handled via simple tie breaking strategies. Argument \a neighborhood
413  specifies the connectivity between points and can be <tt>DirectNeighborhood</tt> (meaning
414  4-neighborhood in 2D and 6-neighborhood in 3D, default) or <tt>IndirectNeighborhood</tt>
415  (meaning 8-neighborhood in 2D and 26-neighborhood in 3D).
416 
417  The watershed variant to be applied can be selected in the \ref WatershedOptions
418  object: When you call <tt>WatershedOptions::regionGrowing()</tt> (default), the flooding
419  algorithm from [1] is used. Alternatively, <tt>WatershedOptions::unionFind()</tt> uses
420  the scan-line algorithm 4.7 from [2]. The latter is faster, but does not support any options
421  (if you pass options nonetheless, they are silently ignored).
422 
423  The region growing algorithm needs a seed for each region. Seeds can either be provided in
424  the destination array \a labels (which will then be overwritten with the result) or computed
425  automatically by an internal call to generateWatershedSeeds(). In the former case you have
426  full control over seed placement, while the latter is more convenient. Automatic seed
427  computation is performed when you provide seeding options via <tt>WatershedOptions::seedOptions()</tt>
428  or when the array \a labels is empty (all zeros), in which case default seeding options
429  are chosen. The destination image should be zero-initialized for automatic seed computation.
430 
431  Further options to be specified via \ref WatershedOptions are:
432 
433  <ul>
434  <li> <tt>keepContours()</tt>: Whether to keep a 1-pixel-wide contour (with label 0) between
435  regions (otherwise, a complete grow is performed, i.e. all pixels are assigned to a region).
436  <li> <tt>stopAtThreshold()</tt>: Whether to stop growing when the boundaryness exceeds a threshold
437  (remaining pixels keep label 0).
438  <li> <tt>biasLabel()</tt>: Whether one region (label) is to be preferred or discouraged by biasing its cost
439  with a given factor (smaller than 1 for preference, larger than 1 for discouragement).
440  </ul>
441 
442  The option <tt>turboAlgorithm()</tt> is implied by method <tt>regionGrowing()</tt> (this is
443  in contrast to watershedsRegionGrowing(), which supports an additional algorithm in 2D only).
444 
445  watershedsMultiArray() returns the number of regions found (= the highest region label, because
446  labels start at 1).
447 
448  <b> Declaration:</b>
449 
450  \code
451  namespace vigra {
452  template <unsigned int N, class T, class S1,
453  class Label, class S2>
454  Label
455  watershedsMultiArray(MultiArrayView<N, T, S1> const & data,
456  MultiArrayView<N, Label, S2> labels, // may also hold input seeds
457  NeighborhoodType neighborhood = DirectNeighborhood,
458  WatershedOptions const & options = WatershedOptions());
459  }
460  \endcode
461 
462  <b> Usage:</b>
463 
464  <b>\#include</b> <vigra/multi_watersheds.hxx><br>
465  Namespace: vigra
466 
467  Example: watersheds of the gradient magnitude (the example works likewise for higher dimensions).
468 
469  \code
470  MultiArray<2, unsigned char> src(Shape2(w, h));
471  ... // read input data
472 
473  // compute gradient magnitude at scale 1.0 as a boundary indicator
474  MultiArray<2, float> gradMag(src.shape());
475  gaussianGradientMagnitude(srcImageRange(src), destImage(gradMag), 1.0);
476 
477  // example 1
478  {
479  // the pixel type of the destination image must be large enough to hold
480  // numbers up to 'max_region_label' to prevent overflow
481  MultiArray<2, unsigned int> labeling(src.shape());
482 
483  // call region-growing algorithm for 4-neighborhood, leave a 1-pixel boundary between
484  // regions, and autogenerate seeds from all gradient minima where the magnitude is
485  // less than 2.0.
486  unsigned int max_region_label =
487  watershedsMultiArray(gradMag, labeling, DirectNeighborhood,
488  WatershedOptions().keepContours()
489  .seedOptions(SeedOptions().minima().threshold(2.0)));
490  }
491 
492  // example 2
493  {
494  MultiArray<2, unsigned int> labeling(src.shape());
495 
496  // compute seeds beforehand (use connected components of all pixels
497  // where the gradient is below 4.0)
498  unsigned int max_region_label = generateWatershedSeeds(gradMag, labeling,
499  SeedOptions().levelSets(4.0));
500 
501  // quantize the gradient image to 256 gray levels
502  float m, M;
503  gradMag.minmax(&m, &M);
504 
505  using namespace multi_math;
506  MultiArray<2, unsigned char> gradMag256(255.0 / (M - m) * (gradMag - m));
507 
508  // call region-growing algorithm with 8-neighborhood,
509  // since the data are 8-bit, a faster priority queue will be used
510  watershedsMultiArray(gradMag256, labeling, IndirectNeighborhood);
511  }
512 
513  // example 3
514  {
515  MultiArray<2, unsigned int> labeling(src.shape());
516 
517  .. // put seeds in 'labeling', e.g. from an interactive labeling program,
518  // make sure that label 1 corresponds to the background
519 
520  // bias the watershed algorithm so that the background is preferred
521  // by reducing the cost for label 1 to 90%
522  watershedsMultiArray(gradMag, labeling,
523  WatershedOptions().biasLabel(1, 0.9));
524  }
525 
526  // example 4
527  {
528  MultiArray<2, unsigned int> labeling(src.shape());
529 
530  // use the fast union-find algorithm with 4-neighborhood
531  watershedsMultiArray(gradMag, labeling, WatershedOptions().unionFind());
532  }
533  \endcode
534 */
536 
537 template <unsigned int N, class T, class S1,
538  class Label, class S2>
539 inline Label
540 watershedsMultiArray(MultiArrayView<N, T, S1> const & data,
541  MultiArrayView<N, Label, S2> labels, // may also hold input seeds
542  NeighborhoodType neighborhood = DirectNeighborhood,
543  WatershedOptions const & options = WatershedOptions())
544 {
545  vigra_precondition(data.shape() == labels.shape(),
546  "watershedsMultiArray(): Shape mismatch between input and output.");
547 
548  GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
549  return lemon_graph::watershedsGraph(graph, data, labels, options);
550 }
551 
552 //@}
553 
554 } // namespace vigra
555 
556 #endif // VIGRA_MULTI_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.10.0 (Thu Jan 8 2015)