36 #ifndef VIGRA_MULTI_WATERSHEDS_HXX
37 #define VIGRA_MULTI_WATERSHEDS_HXX
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"
55 namespace lemon_graph {
57 namespace graph_detail {
59 template <
class Graph,
class T1Map,
class T2Map>
61 prepareWatersheds(Graph
const & g,
63 T2Map & lowestNeighborIndex)
65 typedef typename Graph::NodeIt graph_scanner;
66 typedef typename Graph::OutArcIt neighbor_iterator;
68 for (graph_scanner node(g); node != INVALID; ++node)
70 typename T1Map::value_type lowestValue = data[*node];
71 typename T2Map::value_type lowestIndex = -1;
73 for(neighbor_iterator arc(g, node); arc != INVALID; ++arc)
75 if(data[g.target(*arc)] <= lowestValue)
77 lowestValue = data[g.target(*arc)];
78 lowestIndex = arc.neighborIndex();
81 lowestNeighborIndex[*node] = lowestIndex;
85 template <
class Graph,
class T1Map,
class T2Map,
class T3Map>
86 typename T2Map::value_type
87 unionFindWatersheds(Graph
const & g,
89 T2Map
const & lowestNeighborIndex,
92 typedef typename Graph::NodeIt graph_scanner;
93 typedef typename Graph::OutBackArcIt neighbor_iterator;
94 typedef typename T3Map::value_type LabelType;
96 vigra::detail::UnionFindArray<LabelType> regions;
99 for (graph_scanner node(g); node != INVALID; ++node)
102 LabelType currentLabel = regions.nextFreeLabel();
103 bool hasPlateauNeighbor =
false;
105 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
108 if(lowestNeighborIndex[*node] == arc.neighborIndex() ||
109 lowestNeighborIndex[g.target(*arc)] == g.oppositeIndex(arc.neighborIndex()))
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);
118 if(hasPlateauNeighbor)
121 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
123 if(data[*node] == data[g.target(*arc)])
125 LabelType neighborLabel = regions[labels[g.target(*arc)]];
126 currentLabel = regions.makeUnion(neighborLabel, currentLabel);
132 labels[*node] = regions.finalizeLabel(currentLabel);
135 LabelType count = regions.makeContiguous();
138 for (graph_scanner node(g); node != INVALID; ++node)
140 labels[*node] = regions[labels[*node]];
145 template <
class Graph,
class T1Map,
class T2Map>
146 typename T2Map::value_type
147 generateWatershedSeeds(Graph
const & g,
150 SeedOptions
const & options = SeedOptions())
152 typedef typename T1Map::value_type DataType;
153 typedef unsigned char MarkerType;
155 typename Graph::template NodeMap<MarkerType> minima(g);
157 if(options.mini == SeedOptions::LevelSets)
159 vigra_precondition(options.thresholdIsValid<DataType>(),
160 "generateWatershedSeeds(): SeedOptions.levelSets() must be specified with threshold.");
162 using namespace multi_math;
163 minima = data <= DataType(options.thresh);
167 DataType threshold = options.thresholdIsValid<DataType>()
169 : NumericTraits<DataType>::max();
171 if(options.mini == SeedOptions::ExtendedMinima)
172 extendedLocalMinMaxGraph(g, data, minima, MarkerType(1), threshold,
173 std::less<DataType>(), std::equal_to<DataType>(),
true);
175 localMinMaxGraph(g, data, minima, MarkerType(1), threshold,
176 std::less<DataType>(),
true);
178 return labelGraphWithBackground(g, minima, seeds, MarkerType(0), std::equal_to<MarkerType>());
182 template <
class Graph,
class T1Map,
class T2Map>
183 typename T2Map::value_type
184 seededWatersheds(Graph
const & g,
187 WatershedOptions
const & options)
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;
195 PriorityQueue<Node, CostType, true> pqueue;
197 bool keepContours = ((options.terminate & KeepContours) != 0);
198 LabelType maxRegionLabel = 0;
200 for (graph_scanner node(g); node != INVALID; ++node)
202 LabelType label = labels[*node];
205 if(maxRegionLabel < label)
206 maxRegionLabel = label;
208 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
210 if(labels[g.target(*arc)] == 0)
213 if(label == options.biased_label)
214 pqueue.push(*node, data[*node] * options.bias);
216 pqueue.push(*node, data[*node]);
223 LabelType contourLabel = maxRegionLabel + 1;
226 while(!pqueue.empty())
228 Node node = pqueue.top();
229 CostType cost = pqueue.topPriority();
232 if((options.terminate & StopAtThreshold) && (cost > options.max_cost))
235 LabelType label = labels[node];
237 if(label == contourLabel)
241 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
243 LabelType neighborLabel = labels[g.target(*arc)];
244 if(neighborLabel == 0)
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)];
252 pqueue.push(g.target(*arc), priority);
254 else if(keepContours && (label != neighborLabel) && (neighborLabel != contourLabel))
258 CostType priority = (neighborLabel == options.biased_label)
259 ? data[g.target(*arc)] * options.bias
260 : data[g.target(*arc)];
262 labels[g.target(*arc)] = contourLabel;
270 typename T2Map::iterator k = labels.begin(),
273 if(*k == contourLabel)
277 return maxRegionLabel;
282 template <
class Graph,
class T1Map,
class T2Map>
283 typename T2Map::value_type
284 watershedsGraph(Graph
const & g,
287 WatershedOptions
const & options)
289 if(options.method == WatershedOptions::UnionFind)
291 vigra_precondition(g.maxDegree() <= NumericTraits<unsigned short>::max(),
292 "watershedsGraph(): cannot handle nodes with degree > 65535.");
294 typename Graph::template NodeMap<unsigned short> lowestNeighborIndex(g);
296 graph_detail::prepareWatersheds(g, data, lowestNeighborIndex);
297 return graph_detail::unionFindWatersheds(g, data, lowestNeighborIndex, labels);
299 else if(options.method == WatershedOptions::RegionGrowing)
301 SeedOptions seed_options;
304 if(options.seed_options.mini != SeedOptions::Unspecified)
306 seed_options = options.seed_options;
312 seed_options.mini = SeedOptions::Unspecified;
315 if(seed_options.mini != SeedOptions::Unspecified)
317 graph_detail::generateWatershedSeeds(g, data, labels, seed_options);
320 return graph_detail::seededWatersheds(g, data, labels, options);
324 vigra_precondition(
false,
325 "watershedsGraph(): invalid method in watershed options.");
334 template <
unsigned int N,
class T,
class S1,
335 class Label,
class S2>
337 generateWatershedSeeds(MultiArrayView<N, T, S1>
const & data,
338 MultiArrayView<N, Label, S2> seeds,
340 SeedOptions
const & options = SeedOptions())
342 vigra_precondition(data.shape() == seeds.shape(),
343 "generateWatershedSeeds(): Shape mismatch between input and output.");
345 GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
346 return lemon_graph::graph_detail::generateWatershedSeeds(graph, data, seeds, options);
489 template <
unsigned int N,
class T,
class S1,
490 class Label,
class S2>
493 MultiArrayView<N, Label, S2> labels,
495 WatershedOptions
const & options = WatershedOptions())
497 vigra_precondition(data.shape() == labels.shape(),
498 "watershedsMultiArray(): Shape mismatch between input and output.");
500 GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
501 return lemon_graph::watershedsGraph(graph, data, labels, options);
508 #endif // VIGRA_MULTI_WATERSHEDS_HXX