36 #ifndef VIGRA_MULTI_WATERSHEDS_HXX
37 #define VIGRA_MULTI_WATERSHEDS_HXX
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"
56 namespace lemon_graph {
58 namespace graph_detail {
62 template <
class Graph>
63 struct NeighborIndexFunctor
65 typedef typename Graph::index_type index_type;
67 template <
class NodeIter,
class ArcIter>
68 static index_type
get(Graph
const & g, NodeIter
const &, ArcIter
const & a)
70 return g.id(g.target(*a));
73 template <
class NodeIter,
class ArcIter>
74 static index_type getOpposite(Graph
const & g, NodeIter
const & n, ArcIter
const &)
79 static index_type invalidIndex(Graph
const &)
81 return std::numeric_limits<index_type>::max();
87 template<
unsigned int N,
class DirectedTag>
88 struct NeighborIndexFunctor<GridGraph<N, DirectedTag> >
90 typedef GridGraph<N, DirectedTag> Graph;
93 template <
class NodeIter,
class ArcIter>
94 static index_type
get(Graph
const &, NodeIter
const &, ArcIter
const & a)
96 return a.neighborIndex();
99 template <
class NodeIter,
class ArcIter>
100 static index_type getOpposite(Graph
const & g, NodeIter
const &, ArcIter
const & a)
102 return g.oppositeIndex(a.neighborIndex());
104 static index_type invalidIndex(Graph
const &)
106 return std::numeric_limits<index_type>::max();
110 template <
class Graph,
class T1Map,
class T2Map>
112 prepareWatersheds(Graph
const & g,
114 T2Map & lowestNeighborIndex)
116 typedef typename Graph::NodeIt graph_scanner;
117 typedef typename Graph::OutArcIt neighbor_iterator;
118 typedef NeighborIndexFunctor<Graph> IndexFunctor;
120 for (graph_scanner node(g); node != INVALID; ++node)
122 typename T1Map::value_type lowestValue = data[*node];
123 typename T2Map::value_type lowestIndex = IndexFunctor::invalidIndex(g);
125 for(neighbor_iterator arc(g, *node); arc != INVALID; ++arc)
127 if(data[g.target(*arc)] < lowestValue)
129 lowestValue = data[g.target(*arc)];
133 lowestNeighborIndex[*node] = lowestIndex;
138 template <
class Graph,
class T1Map,
class T2Map,
class T3Map>
139 typename T2Map::value_type
140 unionFindWatersheds(Graph
const & g,
142 T2Map
const & lowestNeighborIndex,
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;
150 vigra::UnionFindArray<LabelType> regions;
153 for (graph_scanner node(g); node != INVALID; ++node)
156 LabelType currentIndex = regions.nextFreeIndex();
158 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
161 if((lowestNeighborIndex[*node] == IndexFunctor::invalidIndex(g) &&
162 lowestNeighborIndex[g.target(*arc)] == IndexFunctor::invalidIndex(g)) ||
164 (lowestNeighborIndex[g.target(*arc)] == IndexFunctor::getOpposite(g, node, arc)))
166 currentIndex = regions.makeUnion(labels[g.target(*arc)], currentIndex);
171 labels[*node] = regions.finalizeIndex(currentIndex);
174 LabelType count = regions.makeContiguous();
177 for (graph_scanner node(g); node != INVALID; ++node)
179 labels[*node] = regions.findLabel(labels[*node]);
184 template <
class Graph,
class T1Map,
class T2Map>
185 typename T2Map::value_type
186 generateWatershedSeeds(Graph
const & g,
189 SeedOptions
const & options = SeedOptions())
191 typedef typename T1Map::value_type DataType;
192 typedef unsigned char MarkerType;
194 typename Graph::template NodeMap<MarkerType> minima(g);
196 if(options.mini == SeedOptions::LevelSets)
198 vigra_precondition(options.thresholdIsValid<DataType>(),
199 "generateWatershedSeeds(): SeedOptions.levelSets() must be specified with threshold.");
201 using namespace multi_math;
202 for(
typename Graph::NodeIt iter(g);iter!=lemon::INVALID;++iter){
203 minima[*iter]= data[*iter] <= DataType(options.thresh);
208 DataType threshold = options.thresholdIsValid<DataType>()
210 : NumericTraits<DataType>::max();
212 if(options.mini == SeedOptions::ExtendedMinima)
213 extendedLocalMinMaxGraph(g, data, minima, MarkerType(1), threshold,
214 std::less<DataType>(), std::equal_to<DataType>(),
true);
216 lemon_graph::localMinMaxGraph(g, data, minima, MarkerType(1), threshold,
217 std::less<DataType>(),
true);
219 return labelGraphWithBackground(g, minima, seeds, MarkerType(0), std::equal_to<MarkerType>());
223 #pragma GCC diagnostic push
224 #pragma GCC diagnostic ignored "-Wsign-compare"
227 template <
class Graph,
class T1Map,
class T2Map>
228 typename T2Map::value_type
229 seededWatersheds(Graph
const & g,
232 WatershedOptions
const & options)
234 typedef typename Graph::Node Node;
235 typedef typename Graph::NodeIt graph_scanner;
236 typedef typename Graph::OutArcIt neighbor_iterator;
237 typedef typename T1Map::value_type CostType;
238 typedef typename T2Map::value_type LabelType;
240 PriorityQueue<Node, CostType, true> pqueue;
242 bool keepContours = ((options.terminate & KeepContours) != 0);
243 LabelType maxRegionLabel = 0;
245 for (graph_scanner node(g); node != INVALID; ++node)
247 LabelType label = labels[*node];
250 if(maxRegionLabel < label)
251 maxRegionLabel = label;
253 for (neighbor_iterator arc(g, *node); arc != INVALID; ++arc)
255 if(labels[g.target(*arc)] == 0)
258 if(label == options.biased_label)
259 pqueue.push(*node, data[*node] * options.bias);
261 pqueue.push(*node, data[*node]);
268 LabelType contourLabel = maxRegionLabel + 1;
271 while(!pqueue.empty())
273 Node node = pqueue.top();
274 CostType cost = pqueue.topPriority();
277 if((options.terminate & StopAtThreshold) && (cost > options.max_cost))
280 LabelType label = labels[node];
282 if(label == contourLabel)
286 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
288 LabelType neighborLabel = labels[g.target(*arc)];
289 if(neighborLabel == 0)
291 labels[g.target(*arc)] = label;
292 CostType priority = (label == options.biased_label)
293 ? data[g.target(*arc)] * options.bias
294 : data[g.target(*arc)];
297 pqueue.push(g.target(*arc), priority);
299 else if(keepContours && (label != neighborLabel) && (neighborLabel != contourLabel))
303 CostType priority = (neighborLabel == options.biased_label)
304 ? data[g.target(*arc)] * options.bias
305 : data[g.target(*arc)];
307 labels[g.target(*arc)] = contourLabel;
321 for(
typename Graph::NodeIt iter(g);iter!=lemon::INVALID;++iter){
322 if(labels[*iter]==contourLabel)
327 return maxRegionLabel;
331 #pragma GCC diagnostic pop
336 template <
class Graph,
class T1Map,
class T2Map>
337 typename T2Map::value_type
338 watershedsGraph(Graph
const & g,
341 WatershedOptions
const & options)
343 if(options.method == WatershedOptions::UnionFind)
345 typedef typename graph_detail::NeighborIndexFunctor<Graph>::index_type index_type;
347 vigra_precondition((index_type)g.maxDegree() <= NumericTraits<index_type>::max(),
348 "watershedsGraph(): cannot handle nodes with degree > 65535.");
350 typename Graph::template NodeMap<index_type> lowestNeighborIndex(g);
352 graph_detail::prepareWatersheds(g, data, lowestNeighborIndex);
353 return graph_detail::unionFindWatersheds(g, data, lowestNeighborIndex, labels);
355 else if(options.method == WatershedOptions::RegionGrowing)
357 SeedOptions seed_options;
360 if(options.seed_options.mini != SeedOptions::Unspecified)
362 seed_options = options.seed_options;
368 seed_options.mini = SeedOptions::Unspecified;
371 if(seed_options.mini != SeedOptions::Unspecified)
373 graph_detail::generateWatershedSeeds(g, data, labels, seed_options);
376 return graph_detail::seededWatersheds(g, data, labels, options);
380 vigra_precondition(
false,
381 "watershedsGraph(): invalid method in watershed options.");
390 template <
unsigned int N,
class T,
class S1,
391 class Label,
class S2>
393 generateWatershedSeeds(MultiArrayView<N, T, S1>
const & data,
394 MultiArrayView<N, Label, S2> seeds,
396 SeedOptions
const & options = SeedOptions())
398 vigra_precondition(data.shape() == seeds.shape(),
399 "generateWatershedSeeds(): Shape mismatch between input and output.");
401 GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
402 return lemon_graph::graph_detail::generateWatershedSeeds(graph, data, seeds, options);
548 template <
unsigned int N,
class T,
class S1,
549 class Label,
class S2>
552 MultiArrayView<N, Label, S2> labels,
554 WatershedOptions
const & options = WatershedOptions())
556 vigra_precondition(data.shape() == labels.shape(),
557 "watershedsMultiArray(): Shape mismatch between input and output.");
559 GridGraph<N, undirected_tag> graph(data.shape(), neighborhood);
560 return lemon_graph::watershedsGraph(graph, data, labels, options);
567 #endif // VIGRA_MULTI_WATERSHEDS_HXX
CoupledHandleCast< TARGET_INDEX, Handle >::reference get(Handle &handle)
Definition: multi_handle.hxx:927
Label watershedsMultiArray(...)
Watershed segmentation of an arbitrary-dimensional array.
detail::SelectIntegerType< 16, detail::UnsignedIntTypes >::type UInt16
16-bit unsigned int
Definition: sized_int.hxx:181
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
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