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

multi_watersheds.hxx |

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 */
538 doxygen_overloaded_function(template <...> Label watershedsMultiArray)
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
vigra::get
CoupledHandleCast< TARGET_INDEX, Handle >::reference get(Handle &handle)
Definition: multi_handle.hxx:927
vigra::UInt16
detail::SelectIntegerType< 16, detail::UnsignedIntTypes >::type UInt16
16-bit unsigned int
Definition: sized_int.hxx:181
vigra
Definition: accessor.hxx:43
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::watershedsMultiArray
Label watershedsMultiArray(...)
Watershed segmentation of an arbitrary-dimensional array.
vigra::IndirectNeighborhood
use direct and indirect neighbors
Definition: multi_fwd.hxx:188
vigra::DirectNeighborhood
use only direct neighbors
Definition: multi_fwd.hxx:187
vigra::NeighborhoodType
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.11.0