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

skeleton.hxx VIGRA

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

© 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