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

skeleton.hxx
1/************************************************************************/
2/* */
3/* Copyright 2013-2014 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_SKELETON_HXX
37#define VIGRA_SKELETON_HXX
38
39#include <vector>
40#include <set>
41#include <map>
42#include "vector_distance.hxx"
43#include "iteratorfacade.hxx"
44#include "pixelneighborhood.hxx"
45#include "graph_algorithms.hxx"
46#include "polygon.hxx"
47
48 namespace vigra
49{
50
51 namespace detail {
52
53 template <class Node>
54 struct SkeletonNode
55{
56 Node parent, principal_child;
57 double length, salience;
58 MultiArrayIndex partial_area;
59 bool is_loop;
60
61 SkeletonNode()
62 : parent(lemon::INVALID)
63 , principal_child(lemon::INVALID)
64 , length(0.0)
65 , salience(1.0)
66 , partial_area(0)
67 , is_loop(false)
68 {}
69
70 SkeletonNode(Node const & s)
71 : parent(s)
72 , principal_child(lemon::INVALID)
73 , length(0.0)
74 , salience(1.0)
75 , partial_area(0)
76 , is_loop(false)
77 {}
78};
79
80 template <class Node>
81 struct SkeletonRegion
82{
83 typedef SkeletonNode<Node> SNode;
84 typedef std::map<Node, SNode> Skeleton;
85
86 Node anchor, lower, upper;
87 Skeleton skeleton;
88
89 SkeletonRegion()
90 : anchor(lemon::INVALID)
91 , lower(NumericTraits<MultiArrayIndex>::max())
92 , upper(NumericTraits<MultiArrayIndex>::min())
93 {}
94
95 void addNode(Node const & n)
96 {
97 skeleton[n] = SNode(n);
98 anchor = n;
99 lower = min(lower, n);
100 upper = max(upper, n);
101 }
102};
103
104 template <class Graph, class Node, class NodeMap>
105 inline unsigned int
106neighborhoodConfiguration(Graph const & g, Node const & n, NodeMap const & labels)
107{
108 typedef typename Graph::OutArcIt ArcIt;
109 typedef typename NodeMap::value_type LabelType;
110
111 LabelType label = labels[n];
112 unsigned int v = 0;
113 for(ArcIt arc(g, n); arc != lemon::INVALID; ++arc)
114 {
115 v = (v << 1) | (labels[g.target(*arc)] == label ? 1 : 0);
116 }
117 return v;
118}
119
120 template <class Node, class Cost>
121 struct SkeletonSimplePoint
122{
123 Node point;
124 Cost cost;
125
126 SkeletonSimplePoint(Node const & p, Cost c)
127 : point(p), cost(c)
128 {}
129
130 bool operator<(SkeletonSimplePoint const & o) const
131 {
132 return cost < o.cost;
133 }
134
135 bool operator>(SkeletonSimplePoint const & o) const
136 {
137 return cost > o.cost;
138 }
139};
140
141 template <class CostMap, class LabelMap>
142 void
143skeletonThinning(CostMap const & cost, LabelMap & labels,
144 bool preserve_endpoints=true)
145{
146 typedef GridGraph<CostMap::actual_dimension> Graph;
147 typedef typename Graph::Node Node;
148 typedef typename Graph::NodeIt NodeIt;
149 typedef typename Graph::OutArcIt ArcIt;
150
151 Graph g(labels.shape(), IndirectNeighborhood);
152 typedef SkeletonSimplePoint<Node, double> SP;
153 // use std::greater because we need the smallest distances at the top of the queue
154 // (std::priority_queue is a max queue by default)
155 std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue;
156
157 bool isSimpleStrong[256] = {
158 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1,
159 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1,
160 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
161 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0,
162 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1,
163 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
164 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
165 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1,
166 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1,
167 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1,
168 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0,
169 };
170
171 bool isSimplePreserveEndPoints[256] = {
172 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1,
173 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1,
174 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
175 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
176 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,
177 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
178 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
179 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
180 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0,
181 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,
182 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
183 };
184
185 bool * isSimplePoint = preserve_endpoints
186 ? isSimplePreserveEndPoints
187 : isSimpleStrong;
188
189 int max_degree = g.maxDegree();
190 double epsilon = 0.5/labels.size(), offset = 0;
191 for (NodeIt node(g); node != lemon::INVALID; ++node)
192 {
193 Node p = *node;
194 if(g.out_degree(p) == max_degree &&
195 labels[p] > 0 &&
196 isSimplePoint[neighborhoodConfiguration(g, p, labels)])
197 {
198 // add an offset to break ties in a FIFI manner
199 pqueue.push(SP(p, cost[p]+offset));
200 offset += epsilon;
201 }
202 }
203
204 while(pqueue.size())
205 {
206 Node p = pqueue.top().point;
207 pqueue.pop();
208
209 if(labels[p] == 0 ||
210 !isSimplePoint[neighborhoodConfiguration(g, p, labels)])
211 {
212 continue; // point already deleted or no longer simple
213 }
214
215 labels[p] = 0; // delete simple point
216
217 for (ArcIt arc(g, p); arc != lemon::INVALID; ++arc)
218 {
219 Node q = g.target(*arc);
220 if(g.out_degree(q) == max_degree &&
221 labels[q] > 0 &&
222 isSimplePoint[neighborhoodConfiguration(g, q, labels)])
223 {
224 pqueue.push(SP(q, cost[q]+offset));
225 offset += epsilon;
226 }
227 }
228 }
229}
230
231 template <class Label, class Labels>
232 struct CheckForHole
233{
234 Label label;
235 Labels const & labels;
236
237 CheckForHole(Label l, Labels const & ls)
238 : label(l)
239 , labels(ls)
240 {}
241
242 template <class Node>
243 bool operator()(Node const & n) const
244 {
245 return labels[n] == label;
246 }
247};
248
249} // namespace detail
250
251/** \addtogroup DistanceTransform
252*/
253//@{
254
255 /** \brief Option object for \ref skeletonizeImage()
256 */
257 struct SkeletonOptions
258{
259 enum SkeletonMode {
260 DontPrune = 0,
261 Prune = 1,
262 Relative = 2,
263 PreserveTopology = 4,
264 Length = 8,
265 Salience = 16,
266 PruneCenterLine = 32,
267 PruneLength = Length + Prune,
268 PruneLengthRelative = PruneLength + Relative,
269 PruneSalience = Salience + Prune,
270 PruneSalienceRelative = PruneSalience + Relative,
271 PruneTopology = PreserveTopology + Prune
272 };
273
274 SkeletonMode mode;
275 double pruning_threshold;
276
277 /** \brief construct with default settings
278
279 (default: <tt>pruneSalienceRelative(0.2, true)</tt>)
280 */
281 SkeletonOptions()
282 : mode(SkeletonMode(PruneSalienceRelative | PreserveTopology))
283 , pruning_threshold(0.2)
284 {}
285
286 /** \brief return the un-pruned skeletong
287 */
288 SkeletonOptions & dontPrune()
289 {
290 mode = DontPrune;
291 return *this;
292 }
293
294 /** \brief return only the region's center line (i.e. skeleton graph diameter)
295 */
296 SkeletonOptions & pruneCenterLine()
297 {
298 mode = PruneCenterLine;
299 return *this;
300 }
301
302 /** \brief Don't prune and return the length of each skeleton segment.
303 */
304 SkeletonOptions & returnLength()
305 {
306 mode = Length;
307 return *this;
308 }
309
310 /** \brief prune skeleton segments whose length is below the given threshold
311
312 If \a preserve_topology is <tt>true</tt> (default), skeleton loops
313 (i.e. parts enclosing a hole in the region) are preserved even if their
314 length is below the threshold. Otherwise, loops are pruned as well.
315 */
316 SkeletonOptions & pruneLength(double threshold, bool preserve_topology=true)
317 {
318 mode = PruneLength;
319 if(preserve_topology)
320 mode = SkeletonMode(mode | PreserveTopology);
321 pruning_threshold = threshold;
322 return *this;
323 }
324
325 /** \brief prune skeleton segments whose relative length is below the given threshold
326
327 This works like <tt>pruneLength()</tt>, but the threshold is specified as a
328 fraction of the maximum segment length in the skeleton.
329 */
330 SkeletonOptions & pruneLengthRelative(double threshold, bool preserve_topology=true)
331 {
332 mode = PruneLengthRelative;
333 if(preserve_topology)
334 mode = SkeletonMode(mode | PreserveTopology);
335 pruning_threshold = threshold;
336 return *this;
337 }
338
339 /** \brief Don't prune and return the salience of each skeleton segment.
340 */
341 SkeletonOptions & returnSalience()
342 {
343 mode = Salience;
344 return *this;
345 }
346
347 /** \brief prune skeleton segments whose salience is below the given threshold
348
349 If \a preserve_topology is <tt>true</tt> (default), skeleton loops
350 (i.e. parts enclosing a hole in the region) are preserved even if their
351 salience is below the threshold. Otherwise, loops are pruned as well.
352 */
353 SkeletonOptions & pruneSalience(double threshold, bool preserve_topology=true)
354 {
355 mode = PruneSalience;
356 if(preserve_topology)
357 mode = SkeletonMode(mode | PreserveTopology);
358 pruning_threshold = threshold;
359 return *this;
360 }
361
362 /** \brief prune skeleton segments whose relative salience is below the given threshold
363
364 This works like <tt>pruneSalience()</tt>, but the threshold is specified as a
365 fraction of the maximum segment salience in the skeleton.
366 */
367 SkeletonOptions & pruneSalienceRelative(double threshold, bool preserve_topology=true)
368 {
369 mode = PruneSalienceRelative;
370 if(preserve_topology)
371 mode = SkeletonMode(mode | PreserveTopology);
372 pruning_threshold = threshold;
373 return *this;
374 }
375
376 /** \brief prune such that only the topology is preserved
377
378 If \a preserve_center is <tt>true</tt> (default), the eccentricity center
379 of the skeleton will not be pruned, even if it is not essential for the topology.
380 Otherwise, the center is only preserved if it is essential. The center is always
381 preserved (and is the only remaining point) when the region has no holes.
382 */
383 SkeletonOptions & pruneTopology(bool preserve_center=true)
384 {
385 if(preserve_center)
386 mode = PruneTopology;
387 else
388 mode = Prune;
389 return *this;
390 }
391};
392
393 template <class T1, class S1,
394 class T2, class S2,
395 class ArrayLike>
396 void
397skeletonizeImageImpl(MultiArrayView<2, T1, S1> const & labels,
398 MultiArrayView<2, T2, S2> dest,
399 ArrayLike * features,
400 SkeletonOptions const & options)
401{
402 static const unsigned int N = 2;
403 typedef typename MultiArrayShape<N>::type Shape;
404 typedef GridGraph<N> Graph;
405 typedef typename Graph::Node Node;
406 typedef typename Graph::NodeIt NodeIt;
407 typedef typename Graph::EdgeIt EdgeIt;
408 typedef typename Graph::OutArcIt ArcIt;
409 typedef typename Graph::OutBackArcIt BackArcIt;
410 typedef double WeightType;
411 typedef detail::SkeletonNode<Node> SNode;
412 typedef std::map<Node, SNode> Skeleton;
413
414 vigra_precondition(labels.shape() == dest.shape(),
415 "skeleton(): shape mismatch between input and output.");
416
417 MultiArray<N, MultiArrayIndex> squared_distance;
418 dest = 0;
419 T1 maxLabel = 0;
420 // find skeleton points
421 {
422 using namespace multi_math;
423
424 MultiArray<N, Shape> vectors(labels.shape());
425 boundaryVectorDistance(labels, vectors, false, OuterBoundary);
426 squared_distance = squaredNorm(vectors);
427
428 ArrayVector<Node> ends_to_be_checked;
429 Graph g(labels.shape());
430 for (EdgeIt edge(g); edge != lemon::INVALID; ++edge)
431 {
432 Node p1 = g.u(*edge),
433 p2 = g.v(*edge);
434 T1 l1 = labels[p1],
435 l2 = labels[p2];
436 maxLabel = max(maxLabel, max(l1, l2));
437 if(l1 == l2)
438 {
439 if(l1 <= 0) // only consider positive labels
440 continue;
441
442 const Node v1 = vectors[p1],
443 v2 = vectors[p2],
444 dp = p2 - p1,
445 dv = v2 - v1 + dp;
446 if(max(abs(dv)) <= 1) // points whose support points coincide or are adjacent
447 continue; // don't belong to the skeleton
448
449 // among p1 and p2, the point which is closer to the bisector
450 // of the support points p1 + v1 and p2 + v2 belongs to the skeleton
451 const MultiArrayIndex d1 = dot(dv, dp),
452 d2 = dot(dv, v1+v2);
453 if(d1*d2 <= 0)
454 {
455 dest[p1] = l1;
456 if(squared_distance[p1] == 4)
457 ends_to_be_checked.push_back(p1);
458 }
459 else
460 {
461 dest[p2] = l2;
462 if(squared_distance[p2] == 4)
463 ends_to_be_checked.push_back(p2);
464 }
465 }
466 else
467 {
468 if(l1 > 0 && max(abs(vectors[p1] + p1 - p2)) > 1)
469 dest[p1] = l1;
470 if(l2 > 0 && max(abs(vectors[p2] + p2 - p1)) > 1)
471 dest[p2] = l2;
472 }
473 }
474
475
476 // add a point when a skeleton line stops short of the shape boundary
477 // FIXME: can this be solved during the initial detection above?
478 Graph g8(labels.shape(), IndirectNeighborhood);
479 for (unsigned k=0; k<ends_to_be_checked.size(); ++k)
480 {
481 // The phenomenon only occurs at points whose distance from the background is 2.
482 // We've put these points into ends_to_be_checked.
483 Node p1 = ends_to_be_checked[k];
484 T2 label = dest[p1];
485 int count = 0;
486 for (ArcIt arc(g8, p1); arc != lemon::INVALID; ++arc)
487 {
488 Node p2 = g8.target(*arc);
489 if(dest[p2] == label && squared_distance[p2] < 4)
490 ++count;
491 }
492 if(count == 0) // point p1 has no neighbor at the boundary => activate one
493 dest[p1+vectors[p1]/2] = label;
494 }
495
496 // from here on, we only need the squared DT, not the vector DT
497 }
498
499 // The true skeleton is defined by the interpixel edges between the
500 // Voronoi regions of the DT. Our skeleton detection algorithm affectively
501 // rounds the interpixel edges to the nearest pixel such that the result
502 // is mainly 8-connected and thin. However, thick skeleton pieces may still
503 // arise when two interpixel contours are only one pixel apart, i.e. a
504 // Voronoi region is only one pixel wide. Since this happens rarely, we
505 // can simply remove these cases by thinning.
506 detail::skeletonThinning(squared_distance, dest);
507
508 if(options.mode == SkeletonOptions::PruneCenterLine)
509 dest = 0;
510
511 // Reduce the full grid graph to a skeleton graph by inserting infinite
512 // edge weights between skeleton pixels and non-skeleton pixels.
513 if(features)
514 features->resize((size_t)maxLabel + 1);
515 ArrayVector<detail::SkeletonRegion<Node> > regions((size_t)maxLabel + 1);
516 Graph g(labels.shape(), IndirectNeighborhood);
517 WeightType maxWeight = g.edgeNum()*sqrt(N),
518 infiniteWeight = 0.5*NumericTraits<WeightType>::max();
519 typename Graph::template EdgeMap<WeightType> weights(g);
520 for (NodeIt node(g); node != lemon::INVALID; ++node)
521 {
522 Node p1 = *node;
523 T2 label = dest[p1];
524 if(label <= 0)
525 continue;
526
527 // FIXME: consider using an AdjacencyListGraph from here on
528 regions[(size_t)label].addNode(p1);
529
530 for (ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
531 {
532 Node p2 = g.target(*arc);
533 if(dest[p2] == label)
534 weights[*arc] = norm(p1-p2);
535 else
536 weights[*arc] = infiniteWeight;
537 }
538 }
539
540 ShortestPathDijkstra<Graph, WeightType> pathFinder(g);
541 // Handle the skeleton of each region individually.
542 for(std::size_t label=1; label < regions.size(); ++label)
543 {
544 Skeleton & skeleton = regions[label].skeleton;
545 if(skeleton.size() == 0) // label doesn't exist
546 continue;
547
548 // Find a diameter (longest path) in the skeleton.
549 Node anchor = regions[label].anchor,
550 lower = regions[label].lower,
551 upper = regions[label].upper + Shape(1);
552
553 pathFinder.run(lower, upper, weights, anchor, lemon::INVALID, maxWeight);
554 anchor = pathFinder.target();
555 pathFinder.reRun(weights, anchor, lemon::INVALID, maxWeight);
556 anchor = pathFinder.target();
557
558 Polygon<Shape> center_line;
559 center_line.push_back_unsafe(anchor);
560 while(pathFinder.predecessors()[center_line.back()] != center_line.back())
561 center_line.push_back_unsafe(pathFinder.predecessors()[center_line.back()]);
562
563 if(options.mode == SkeletonOptions::PruneCenterLine)
564 {
565 for(unsigned int k=0; k<center_line.size(); ++k)
566 dest[center_line[k]] = (T2)label;
567 continue; // to next label
568 }
569
570 // Perform the eccentricity transform of the skeleton
571 Node center = center_line[roundi(center_line.arcLengthQuantile(0.5))];
572 pathFinder.reRun(weights, center, lemon::INVALID, maxWeight);
573
574 bool compute_salience = (options.mode & SkeletonOptions::Salience) != 0;
575 ArrayVector<Node> raw_skeleton(pathFinder.discoveryOrder());
576 // from periphery to center: create skeleton tree and compute salience
577 for(int k=raw_skeleton.size()-1; k >= 0; --k)
578 {
579 Node p1 = raw_skeleton[k],
580 p2 = pathFinder.predecessors()[p1];
581 SNode & n1 = skeleton[p1];
582 SNode & n2 = skeleton[p2];
583 n1.parent = p2;
584
585
586 // remove non-skeleton edges (i.e. set weight = infiniteWeight)
587 for (BackArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
588 {
589 Node p = g.target(*arc);
590 if(weights[*arc] == infiniteWeight)
591 continue; // edge never was in the graph
592 if(p == p2)
593 continue; // edge belongs to the skeleton
594 if(pathFinder.predecessors()[p] == p1)
595 continue; // edge belongs to the skeleton
596 if(n1.principal_child == lemon::INVALID ||
597 skeleton[p].principal_child == lemon::INVALID)
598 continue; // edge may belong to a loop => test later
599 weights[*arc] = infiniteWeight;
600 }
601
602 // propagate length to parent if this is the longest subtree
603 WeightType l = n1.length + norm(p1-p2);
604 if(n2.length < l)
605 {
606 n2.length = l;
607 n2.principal_child = p1;
608 }
609
610 if(compute_salience)
611 {
612 const double min_length = 4.0; // salience is meaningless for shorter segments due
613 // to quantization noise (staircasing) of the boundary
614 if(n1.length >= min_length)
615 {
616 n1.salience = max(n1.salience, (n1.length + 0.5) / sqrt(squared_distance[p1]));
617
618 // propagate salience to parent if this is the most salient subtree
619 if(n2.salience < n1.salience)
620 n2.salience = n1.salience;
621 }
622 }
623 else if(options.mode == SkeletonOptions::DontPrune)
624 n1.salience = dest[p1];
625 else
626 n1.salience = n1.length;
627 }
628
629 // from center to periphery: propagate salience and compute twice the partial area
630 for(int k=0; k < (int)raw_skeleton.size(); ++k)
631 {
632 Node p1 = raw_skeleton[k];
633 SNode & n1 = skeleton[p1];
634 Node p2 = n1.parent;
635 SNode & n2 = skeleton[p2];
636
637 if(p1 == n2.principal_child)
638 {
639 n1.length = n2.length;
640 n1.salience = n2.salience;
641 }
642 else
643 {
644 n1.length += norm(p2-p1);
645 }
646 n1.partial_area = n2.partial_area + (p1[0]*p2[1] - p1[1]*p2[0]);
647 }
648
649 // always treat eccentricity center as a loop, so that it cannot be pruned
650 // away unless (options.mode & PreserveTopology) is false.
651 skeleton[center].is_loop = true;
652
653 // from periphery to center: * find and propagate loops
654 // * delete branches not reaching the boundary
655 detail::CheckForHole<std::size_t, MultiArrayView<2, T1, S1> > hasNoHole(label, labels);
656 int hole_count = 0;
657 double total_length = 0.0;
658 for(int k=raw_skeleton.size()-1; k >= 0; --k)
659 {
660 Node p1 = raw_skeleton[k];
661 SNode & n1 = skeleton[p1];
662
663 if(n1.principal_child == lemon::INVALID)
664 {
665 for (ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
666 {
667 Node p2 = g.target(*arc);
668 SNode * n2 = &skeleton[p2];
669
670 if(n1.parent == p2)
671 continue; // going back to the parent can't result in a loop
672 if(weights[*arc] == infiniteWeight)
673 continue; // p2 is not in the tree or the loop has already been handled
674 // compute twice the area exclosed by the potential loop
675 MultiArrayIndex area2 = abs(n1.partial_area - (p1[0]*p2[1] - p1[1]*p2[0]) - n2->partial_area);
676 if(area2 <= 3) // area is too small to enclose a hole => loop is a discretization artifact
677 continue;
678
679 // use Dijkstra to find the loop
680 weights[*arc] = infiniteWeight;
681 pathFinder.reRun(weights, p1, p2);
682 Polygon<Shape2> poly;
683 {
684 poly.push_back_unsafe(p1);
685 poly.push_back_unsafe(p2);
686 Node p = p2;
687 do
688 {
689 p = pathFinder.predecessors()[p];
690 poly.push_back_unsafe(p);
691 }
692 while(p != pathFinder.predecessors()[p]);
693 }
694 // check if the loop contains a hole or is just a discretization artifact
695 if(!inspectPolygon(poly, hasNoHole))
696 {
697 // it's a genuine loop => mark it and propagate salience
698 ++hole_count;
699 total_length += n1.length + n2->length;
700 double max_salience = max(n1.salience, n2->salience);
701 for(decltype(poly.size()) p=1; p<poly.size(); ++p)
702 {
703 SNode & n = skeleton[poly[p]];
704 n.is_loop = true;
705 n.salience = max(n.salience, max_salience);
706 }
707 }
708 }
709 // delete skeleton branches that are not loops and don't reach the shape border
710 // (these branches are discretization artifacts)
711 if(!n1.is_loop && squared_distance[p1] >= 4)
712 {
713 SNode * n = &n1;
714 while(true)
715 {
716 n->salience = 0;
717 // remove all of p1's edges
718 for(ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
719 {
720 weights[*arc] = infiniteWeight;
721 }
722 if(skeleton[n->parent].principal_child != p1)
723 break;
724 p1 = n->parent;
725 n = &skeleton[p1];
726 }
727 }
728 }
729
730 if(n1.is_loop)
731 skeleton[n1.parent].is_loop = true;
732 }
733
734 bool dont_prune = (options.mode & SkeletonOptions::Prune) == 0;
735 bool preserve_topology = (options.mode & SkeletonOptions::PreserveTopology) != 0 ||
736 options.mode == SkeletonOptions::Prune;
737 bool relative_pruning = (options.mode & SkeletonOptions::Relative) != 0;
738 WeightType threshold = (options.mode == SkeletonOptions::PruneTopology ||
739 options.mode == SkeletonOptions::Prune)
740 ? infiniteWeight
741 : relative_pruning
742 ? options.pruning_threshold*skeleton[center].salience
743 : options.pruning_threshold;
744 // from center to periphery: write result
745 int branch_count = 0;
746 double average_length = 0;
747 for(int k=0; k < (int)raw_skeleton.size(); ++k)
748 {
749 Node p1 = raw_skeleton[k];
750 SNode & n1 = skeleton[p1];
751 if(n1.principal_child == lemon::INVALID &&
752 n1.salience >= threshold &&
753 !n1.is_loop)
754 {
755 ++branch_count;
756 average_length += n1.length;
757 total_length += n1.length;
758 }
759 if(dont_prune)
760 dest[p1] = n1.salience;
761 else if(preserve_topology)
762 {
763 if(!n1.is_loop && n1.salience < threshold)
764 dest[p1] = 0;
765 }
766 else if(p1 != center && n1.salience < threshold)
767 dest[p1] = 0;
768 }
769 if(branch_count > 0)
770 average_length /= branch_count;
771
772 if(features)
773 {
774 (*features)[label].diameter = center_line.length();
775 (*features)[label].total_length = total_length;
776 (*features)[label].average_length = average_length;
777 (*features)[label].branch_count = branch_count;
778 (*features)[label].hole_count = hole_count;
779 (*features)[label].center = center;
780 (*features)[label].terminal1 = center_line.front();
781 (*features)[label].terminal2 = center_line.back();
782 (*features)[label].euclidean_diameter = norm(center_line.back()-center_line.front());
783 }
784 }
785
786 if(options.mode == SkeletonOptions::Prune)
787 detail::skeletonThinning(squared_distance, dest, false);
788}
789
790 class SkeletonFeatures
791{
792 public:
793 double diameter, total_length, average_length, euclidean_diameter;
794 UInt32 branch_count, hole_count;
795 Shape2 center, terminal1, terminal2;
796
797 SkeletonFeatures()
798 : diameter(0)
799 , total_length(0)
800 , average_length(0)
801 , euclidean_diameter(0)
802 , branch_count(0)
803 , hole_count(0)
804 {}
805};
806
807/********************************************************/
808/* */
809/* skeletonizeImage */
810/* */
811/********************************************************/
812
813 /*
814 To compute the skeleton reliably in higher dimensions, we have to work on
815 a topological grid. The tricks to work with rounded skeletons on the
816 pixel grid probably don't generalize from 2D to 3D and higher. Specifically:
817
818 * Compute Voronoi regions of the vector distance transformation according to
819 identical support point to make sure that disconnected Voronoi regions
820 still get only a single label.
821 * Merge Voronoi regions whose support points are adjacent.
822 * Mark skeleton candidates on the interpixel grid after the basic merge.
823 * Detect skeleton segments simply by connected components labeling in the interpixel grid.
824 * Skeleton segments form hyperplanes => use this property to compute segment
825 attributes.
826 * Detect holes (and therefore, skeleton segments that are critical for topology)
827 by computing the depth of each region/surface in the homotopy tree.
828 * Add a pruning mode where holes are only preserved if their size exceeds a threshold.
829
830 To implement this cleanly, we first need a good implementation of the topological grid graph.
831 */
832// template <unsigned int N, class T1, class S1,
833 // class T2, class S2>
834// void
835// skeletonizeImage(MultiArrayView<N, T1, S1> const & labels,
836 // MultiArrayView<N, T2, S2> dest,
837 // SkeletonOptions const & options = SkeletonOptions())
838// {
839
840 /** \brief Skeletonization of all regions in a labeled 2D image.
841
842 <b> Declarations:</b>
843
844 \code
845 namespace vigra {
846 template <class T1, class S1,
847 class T2, class S2>
848 void
849 skeletonizeImage(MultiArrayView<2, T1, S1> const & labels,
850 MultiArrayView<2, T2, S2> dest,
851 SkeletonOptions const & options = SkeletonOptions());
852 }
853 \endcode
854
855 This function computes the skeleton for each region in the 2D label image \a labels
856 and paints the results into the result image \a dest. Input label
857 <tt>0</tt> is interpreted as background and always ignored. Skeletons will be
858 marked with the same label as the corresponding region (unless options
859 <tt>returnLength()</tt> or <tt>returnSalience()</tt> are selected, see below).
860 Non-skeleton pixels will receive label <tt>0</tt>.
861
862 For each region, the algorithm proceeds in the following steps:
863 <ol>
864 <li>Compute the \ref boundaryVectorDistance() relative to the \ref OuterBoundary of the region.</li>
865 <li>Mark the raw skeleton: find 4-adjacent pixel pairs whose nearest boundary points are neither equal
866 nor adjacent and mark one pixel of the pair as a skeleton candidate. The resulting raw skeleton
867 is 8-connected and thin. Skip the remaining steps when option <tt>dontPrune()</tt> is selected.</li>
868 <li>Compute the eccentricity transform of the raw skeleton and turn the skeleton into a tree
869 whose root is the eccentricity center. When option <tt>pruneCenterLine()</tt> is selected,
870 delete all skeleton points that do not belong to the two longest tree branches and
871 skip the remaining steps.</li>
872 <li>For each pixel on the skeleton, compute its <tt>length</tt> attribute as the depth of the
873 pixel's longest subtree. Compute its <tt>salience</tt> attribute as the ratio between
874 <tt>length</tt> and <tt>distance</tt>, where <tt>distance</tt> is the pixel's distance to
875 the nearest boundary point according to the distance transform. It holds that <tt>length >= 0.5</tt>
876 and <tt>salience >= 1.0</tt>.</li>
877 <li>Detect skeleton branching points and define <i>skeleton segments</i> as maximal connected pieces
878 without branching points.</li>
879 <li>Compute <tt>length</tt> and <tt>salience</tt> of each segment as the maximum of these
880 attributes among the pixels in the segment. When options <tt>returnLength()</tt> or
881 <tt>returnSalience()</tt> are selected, skip the remaining steps and return the
882 requested segment attribute in <tt>dest</tt>. In this case, <tt>dest</tt>'s
883 <tt>value_type</tt> should be a floating point type to exactly accomodate the
884 attribute values.</li>
885 <li>Detect minimal cycles in the raw skeleton that enclose holes in the region (if any) and mark
886 the corresponding pixels as critical for skeleton topology.</li>
887 <li>Prune skeleton segments according to the selected pruning strategy and return the result.
888 The following pruning strategies are available:
889 <ul>
890 <li><tt>pruneLength(threshold, preserve_topology)</tt>: Retain only segments whose length attribute
891 exceeds the given <tt>threshold</tt>. When <tt>preserve_topology</tt> is true
892 (the defult), cycles around holes are preserved regardless of their length.
893 Otherwise, they are pruned as well.</li>
894 <li><tt>pruneLengthRelative(threshold, preserve_topology)</tt>: Like <tt>pruneLength()</tt>,
895 but the threshold is specified as a fraction of the maximum segment length in
896 the present region.</li>
897 <li><tt>pruneSalience(threshold, preserve_topology)</tt>: Retain only segments whose salience attribute
898 exceeds the given <tt>threshold</tt>. When <tt>preserve_topology</tt> is true
899 (the defult), cycles around holes are preserved regardless of their salience.
900 Otherwise, they are pruned as well.</li>
901 <li><tt>pruneSalienceRelative(threshold, preserve_topology)</tt>: Like <tt>pruneSalience()</tt>,
902 but the threshold is specified as a fraction of the maximum segment salience in
903 the present region.</li>
904 <li><tt>pruneTopology(preserve_center)</tt>: Retain only segments that are essential for the region's
905 topology. If <tt>preserve_center</tt> is true (the default), the eccentricity
906 center is also preserved, even if it is not essential. Otherwise, it might be
907 removed. The eccentricity center is always the only remaining point when
908 the region has no holes.</li>
909 </ul></li>
910 </ol>
911
912 The skeleton has the following properties:
913 <ul>
914 <li>It is 8-connected and thin (except when two independent branches happen to run alongside
915 before they divert). Skeleton points are defined by rounding the exact Euclidean skeleton
916 locations to the nearest pixel.</li>
917 <li>Skeleton branches terminate either at the region boundary or at a cycle. There are no branch
918 end points in the region interior.</li>
919 <li>The salience threshold acts as a scale parameter: Large thresholds only retain skeleton
920 branches characterizing the general region shape. When the threshold gets smaller, ever
921 more detailed boundary bulges will be represented by a skeleton branch.</li>
922 </ul>
923
924 Remark: If you have an application where a skeleton graph would be more useful
925 than a skeleton image, function <tt>skeletonizeImage()</tt> can be changed/extended easily.
926
927 <b> Usage:</b>
928
929 <b>\#include</b> <vigra/skeleton.hxx><br/>
930 Namespace: vigra
931
932 \code
933 Shape2 shape(width, height);
934 MultiArray<2, UInt32> source(shape);
935 MultiArray<2, UInt32> dest(shape);
936 ...
937
938 // Skeletonize and keep only those segments that are at least 10% of the maximum
939 // length (the maximum length is half the skeleton diameter).
940 skeletonizeImage(source, dest,
941 SkeletonOptions().pruneLengthRelative(0.1));
942 \endcode
943
944 \see vigra::boundaryVectorDistance()
945 */
946doxygen_overloaded_function(template <...> void skeletonizeImage)
947
948 template <class T1, class S1,
949 class T2, class S2>
950 void
951 skeletonizeImage(MultiArrayView<2, T1, S1> const & labels,
952 MultiArrayView<2, T2, S2> dest,
953 SkeletonOptions const & options = SkeletonOptions())
954{
955 skeletonizeImageImpl(labels, dest, (ArrayVector<SkeletonFeatures>*)0, options);
956}
957
958 template <class T, class S>
959 void
960extractSkeletonFeatures(MultiArrayView<2, T, S> const & labels,
961 ArrayVector<SkeletonFeatures> & features,
962 SkeletonOptions const & options = SkeletonOptions())
963{
964 MultiArray<2, float> skeleton(labels.shape());
965 skeletonizeImageImpl(labels, skeleton, &features, options);
966}
967
968//@}
969
970} //-- namespace vigra
971
972#endif //-- VIGRA_SKELETON_HXX
vigra::MultiArrayShape::type
TinyVector< MultiArrayIndex, N > type
Definition multi_shape.hxx:272
vigra::RGBValue
Class for a single RGB value.
Definition rgbvalue.hxx:128
vigra
vigra::OuterBoundary
@ OuterBoundary
Pixels just outside of each region.
Definition multi_distance.hxx:836
vigra::UInt32
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition sized_int.hxx:183
vigra::skeletonizeImage
void skeletonizeImage(...)
Skeletonization of all regions in a labeled 2D image.
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
vigra::roundi
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition fixedpoint.hxx:1775
vigra::IndirectNeighborhood
@ IndirectNeighborhood
use direct and indirect neighbors
Definition multi_fwd.hxx:188
vigra::operator>
bool operator>(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
greater
Definition fixedpoint.hxx:530
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
vigra::MultiArrayIndex
std::ptrdiff_t MultiArrayIndex
Definition multi_fwd.hxx:60
vigra::operator<
bool operator<(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
less than
Definition fixedpoint.hxx:512
vigra::boundaryVectorDistance
void boundaryVectorDistance(...)
Compute the vector distance transform to the implicit boundaries of a multi-dimensional label array.
vigra::SkeletonOptions
Option object for skeletonizeImage()
Definition skeleton.hxx:258
vigra::SkeletonOptions::pruneSalience
SkeletonOptions & pruneSalience(double threshold, bool preserve_topology=true)
prune skeleton segments whose salience is below the given threshold
Definition skeleton.hxx:353
vigra::SkeletonOptions::returnSalience
SkeletonOptions & returnSalience()
Don't prune and return the salience of each skeleton segment.
Definition skeleton.hxx:341
vigra::SkeletonOptions::pruneLength
SkeletonOptions & pruneLength(double threshold, bool preserve_topology=true)
prune skeleton segments whose length is below the given threshold
Definition skeleton.hxx:316
vigra::SkeletonOptions::returnLength
SkeletonOptions & returnLength()
Don't prune and return the length of each skeleton segment.
Definition skeleton.hxx:304
vigra::SkeletonOptions::pruneCenterLine
SkeletonOptions & pruneCenterLine()
return only the region's center line (i.e. skeleton graph diameter)
Definition skeleton.hxx:296
vigra::SkeletonOptions::SkeletonOptions
SkeletonOptions()
construct with default settings
Definition skeleton.hxx:281
vigra::SkeletonOptions::pruneSalienceRelative
SkeletonOptions & pruneSalienceRelative(double threshold, bool preserve_topology=true)
prune skeleton segments whose relative salience is below the given threshold
Definition skeleton.hxx:367
vigra::SkeletonOptions::pruneTopology
SkeletonOptions & pruneTopology(bool preserve_center=true)
prune such that only the topology is preserved
Definition skeleton.hxx:383
vigra::SkeletonOptions::pruneLengthRelative
SkeletonOptions & pruneLengthRelative(double threshold, bool preserve_topology=true)
prune skeleton segments whose relative length is below the given threshold
Definition skeleton.hxx:330
vigra::SkeletonOptions::dontPrune
SkeletonOptions & dontPrune()
return the un-pruned skeletong
Definition skeleton.hxx:288

© 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.12.2