36 #ifndef VIGRA_SKELETON_HXX
37 #define VIGRA_SKELETON_HXX
42 #include "vector_distance.hxx"
43 #include "iteratorfacade.hxx"
44 #include "pixelneighborhood.hxx"
45 #include "graph_algorithms.hxx"
55 Node parent, principal_child;
56 double length, salience;
61 : parent(lemon::INVALID)
62 , principal_child(lemon::INVALID)
69 SkeletonNode(Node
const & s)
71 , principal_child(lemon::INVALID)
82 typedef SkeletonNode<Node> SNode;
83 typedef std::map<Node, SNode> Skeleton;
85 Node anchor, lower, upper;
89 : anchor(lemon::INVALID)
94 void addNode(Node
const & n)
96 skeleton[n] = SNode(n);
98 lower = min(lower, n);
99 upper = max(upper, n);
103 template <
class Graph,
class Node,
class NodeMap>
105 neighborhoodConfiguration(Graph
const & g, Node
const & n, NodeMap
const & labels)
107 typedef typename Graph::OutArcIt ArcIt;
108 typedef typename NodeMap::value_type LabelType;
110 LabelType label = labels[n];
112 for(ArcIt arc(g, n); arc != lemon::INVALID; ++arc)
114 v = (v << 1) | (labels[g.target(*arc)] == label ? 1 : 0);
119 template <
class Node,
class Cost>
120 struct SkeletonSimplePoint
125 SkeletonSimplePoint(Node
const & p, Cost c)
129 bool operator<(SkeletonSimplePoint
const & o)
const
131 return cost < o.cost;
134 bool operator>(SkeletonSimplePoint
const & o)
const
136 return cost > o.cost;
140 template <
class CostMap,
class LabelMap>
142 skeletonThinning(CostMap
const & cost, LabelMap & labels,
143 bool preserve_endpoints=
true)
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;
151 typedef SkeletonSimplePoint<Node, double> SP;
154 std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue;
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,
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
184 bool * isSimplePoint = preserve_endpoints
185 ? isSimplePreserveEndPoints
188 int max_degree = g.maxDegree();
189 double epsilon = 0.5/labels.size(), offset = 0;
190 for (NodeIt node(g); node != lemon::INVALID; ++node)
193 if(g.out_degree(p) == max_degree &&
195 isSimplePoint[neighborhoodConfiguration(g, p, labels)])
198 pqueue.push(SP(p, cost[p]+offset));
205 Node p = pqueue.top().point;
209 !isSimplePoint[neighborhoodConfiguration(g, p, labels)])
216 for (ArcIt arc(g, p); arc != lemon::INVALID; ++arc)
218 Node q = g.target(*arc);
219 if(g.out_degree(q) == max_degree &&
221 isSimplePoint[neighborhoodConfiguration(g, q, labels)])
223 pqueue.push(SP(q, cost[q]+offset));
230 template <
class Label,
class Labels>
234 Labels
const & labels;
236 CheckForHole(Label l, Labels
const & ls)
241 template <
class Node>
242 bool operator()(Node
const & n)
const
244 return labels[n] == label;
262 PreserveTopology = 4,
265 PruneCenterLine = 32,
266 PruneLength = Length + Prune,
267 PruneLengthRelative = PruneLength + Relative,
268 PruneSalience = Salience + Prune,
269 PruneSalienceRelative = PruneSalience + Relative,
270 PruneTopology = PreserveTopology + Prune
274 double pruning_threshold;
281 : mode(SkeletonMode(PruneSalienceRelative | PreserveTopology))
282 , pruning_threshold(0.2)
297 mode = PruneCenterLine;
318 if(preserve_topology)
319 mode = SkeletonMode(mode | PreserveTopology);
320 pruning_threshold = threshold;
331 mode = PruneLengthRelative;
332 if(preserve_topology)
333 mode = SkeletonMode(mode | PreserveTopology);
334 pruning_threshold = threshold;
354 mode = PruneSalience;
355 if(preserve_topology)
356 mode = SkeletonMode(mode | PreserveTopology);
357 pruning_threshold = threshold;
368 mode = PruneSalienceRelative;
369 if(preserve_topology)
370 mode = SkeletonMode(mode | PreserveTopology);
371 pruning_threshold = threshold;
385 mode = PruneTopology;
392 template <
class T1,
class S1,
396 skeletonizeImageImpl(MultiArrayView<2, T1, S1>
const & labels,
397 MultiArrayView<2, T2, S2> dest,
398 ArrayLike * features,
399 SkeletonOptions
const & options)
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;
413 vigra_precondition(labels.shape() == dest.shape(),
414 "skeleton(): shape mismatch between input and output.");
416 MultiArray<N, MultiArrayIndex> squared_distance;
421 using namespace multi_math;
423 MultiArray<N, Shape> vectors(labels.shape());
427 ArrayVector<Node> ends_to_be_checked;
428 Graph g(labels.shape());
431 Node p1 = g.u(*
edge),
435 maxLabel = max(maxLabel, max(l1, l2));
441 const Node v1 = vectors[p1],
445 if(max(
abs(dv)) <= 1)
455 if(squared_distance[p1] == 4)
456 ends_to_be_checked.push_back(p1);
461 if(squared_distance[p2] == 4)
462 ends_to_be_checked.push_back(p2);
467 if(l1 > 0 && max(
abs(vectors[p1] + p1 - p2)) > 1)
469 if(l2 > 0 && max(
abs(vectors[p2] + p2 - p1)) > 1)
478 for (
unsigned k=0; k<ends_to_be_checked.size(); ++k)
482 Node p1 = ends_to_be_checked[k];
485 for (ArcIt arc(g8, p1); arc != lemon::INVALID; ++arc)
487 Node p2 = g8.target(*arc);
488 if(dest[p2] == label && squared_distance[p2] < 4)
492 dest[p1+vectors[p1]/2] = label;
505 detail::skeletonThinning(squared_distance, dest);
507 if(options.mode == SkeletonOptions::PruneCenterLine)
513 features->resize((
size_t)maxLabel + 1);
514 ArrayVector<detail::SkeletonRegion<Node> > regions((
size_t)maxLabel + 1);
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)
527 regions[(size_t)label].addNode(p1);
529 for (ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
531 Node p2 = g.target(*arc);
532 if(dest[p2] == label)
533 weights[*arc] =
norm(p1-p2);
535 weights[*arc] = infiniteWeight;
539 ShortestPathDijkstra<Graph, WeightType> pathFinder(g);
541 for(std::size_t label=1; label < regions.size(); ++label)
543 Skeleton & skeleton = regions[label].skeleton;
544 if(skeleton.size() == 0)
548 Node anchor = regions[label].anchor,
549 lower = regions[label].lower,
550 upper = regions[label].upper + Shape(1);
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();
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()]);
562 if(options.mode == SkeletonOptions::PruneCenterLine)
564 for(
unsigned int k=0; k<center_line.size(); ++k)
565 dest[center_line[k]] = (T2)label;
570 Node center = center_line[
roundi(center_line.arcLengthQuantile(0.5))];
571 pathFinder.reRun(weights, center, lemon::INVALID, maxWeight);
573 bool compute_salience = (options.mode & SkeletonOptions::Salience) != 0;
574 ArrayVector<Node> raw_skeleton(pathFinder.discoveryOrder());
576 for(
int k=raw_skeleton.size()-1; k >= 0; --k)
578 Node p1 = raw_skeleton[k],
579 p2 = pathFinder.predecessors()[p1];
580 SNode & n1 = skeleton[p1];
581 SNode & n2 = skeleton[p2];
586 for (BackArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
588 Node p = g.target(*arc);
589 if(weights[*arc] == infiniteWeight)
593 if(pathFinder.predecessors()[p] == p1)
595 if(n1.principal_child == lemon::INVALID ||
596 skeleton[p].principal_child == lemon::INVALID)
598 weights[*arc] = infiniteWeight;
602 WeightType l = n1.length +
norm(p1-p2);
606 n2.principal_child = p1;
611 const double min_length = 4.0;
613 if(n1.length >= min_length)
615 n1.salience = max(n1.salience, (n1.length + 0.5) /
sqrt(squared_distance[p1]));
618 if(n2.salience < n1.salience)
619 n2.salience = n1.salience;
622 else if(options.mode == SkeletonOptions::DontPrune)
623 n1.salience = dest[p1];
625 n1.salience = n1.length;
629 for(
int k=0; k < (int)raw_skeleton.size(); ++k)
631 Node p1 = raw_skeleton[k];
632 SNode & n1 = skeleton[p1];
634 SNode & n2 = skeleton[p2];
636 if(p1 == n2.principal_child)
638 n1.length = n2.length;
639 n1.salience = n2.salience;
643 n1.length +=
norm(p2-p1);
645 n1.partial_area = n2.partial_area + (p1[0]*p2[1] - p1[1]*p2[0]);
650 skeleton[center].is_loop =
true;
654 detail::CheckForHole<std::size_t, MultiArrayView<2, T1, S1> > hasNoHole(label, labels);
656 double total_length = 0.0;
657 for(
int k=raw_skeleton.size()-1; k >= 0; --k)
659 Node p1 = raw_skeleton[k];
660 SNode & n1 = skeleton[p1];
662 if(n1.principal_child == lemon::INVALID)
664 for (ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
666 Node p2 = g.target(*arc);
667 SNode * n2 = &skeleton[p2];
671 if(weights[*arc] == infiniteWeight)
674 MultiArrayIndex area2 =
abs(n1.partial_area - (p1[0]*p2[1] - p1[1]*p2[0]) - n2->partial_area);
679 WeightType edge_length = weights[*arc];
680 weights[*arc] = infiniteWeight;
681 pathFinder.reRun(weights, p1, p2);
682 Polygon<Shape2> poly;
684 poly.push_back_unsafe(p1);
685 poly.push_back_unsafe(p2);
689 p = pathFinder.predecessors()[p];
690 poly.push_back_unsafe(p);
692 while(p != pathFinder.predecessors()[p]);
695 if(!inspectPolygon(poly, hasNoHole))
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)
703 SNode & n = skeleton[poly[p]];
705 n.salience = max(n.salience, max_salience);
711 if(!n1.is_loop && squared_distance[p1] >= 4)
718 for(ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
720 weights[*arc] = infiniteWeight;
722 if(skeleton[n->parent].principal_child != p1)
731 skeleton[n1.parent].is_loop =
true;
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)
742 ? options.pruning_threshold*skeleton[center].salience
743 : options.pruning_threshold;
745 int branch_count = 0;
746 double average_length = 0;
747 for(
int k=0; k < (
int)raw_skeleton.size(); ++k)
749 Node p1 = raw_skeleton[k];
750 SNode & n1 = skeleton[p1];
752 SNode & n2 = skeleton[p2];
753 if(n1.principal_child == lemon::INVALID &&
754 n1.salience >= threshold &&
758 average_length += n1.length;
759 total_length += n1.length;
762 dest[p1] = n1.salience;
763 else if(preserve_topology)
765 if(!n1.is_loop && n1.salience < threshold)
768 else if(p1 != center && n1.salience < threshold)
772 average_length /= branch_count;
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());
788 if(options.mode == SkeletonOptions::Prune)
789 detail::skeletonThinning(squared_distance, dest,
false);
792 class SkeletonFeatures
795 double diameter, total_length, average_length, euclidean_diameter;
796 UInt32 branch_count, hole_count;
797 Shape2 center, terminal1, terminal2;
803 , euclidean_diameter(0)
950 template <
class T1,
class S1,
954 MultiArrayView<2, T2, S2> dest,
955 SkeletonOptions
const & options = SkeletonOptions())
957 skeletonizeImageImpl(labels, dest, (ArrayVector<SkeletonFeatures>*)0, options);
960 template <
class T,
class S>
962 extractSkeletonFeatures(MultiArrayView<2, T, S>
const & labels,
963 ArrayVector<SkeletonFeatures> & features,
964 SkeletonOptions
const & options = SkeletonOptions())
966 MultiArray<2, float> skeleton(labels.shape());
967 skeletonizeImageImpl(labels, skeleton, &features, options);
974 #endif //-- VIGRA_SKELETON_HXX