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

© 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 (Fri Feb 21 2014)