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 weights[*arc] = infiniteWeight;
680 pathFinder.reRun(weights, p1, p2);
681 Polygon<Shape2> poly;
683 poly.push_back_unsafe(p1);
684 poly.push_back_unsafe(p2);
688 p = pathFinder.predecessors()[p];
689 poly.push_back_unsafe(p);
691 while(p != pathFinder.predecessors()[p]);
694 if(!inspectPolygon(poly, hasNoHole))
698 total_length += n1.length + n2->length;
699 double max_salience = max(n1.salience, n2->salience);
700 for(decltype(poly.size()) p=1; p<poly.size(); ++p)
702 SNode & n = skeleton[poly[p]];
704 n.salience = max(n.salience, max_salience);
710 if(!n1.is_loop && squared_distance[p1] >= 4)
717 for(ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
719 weights[*arc] = infiniteWeight;
721 if(skeleton[n->parent].principal_child != p1)
730 skeleton[n1.parent].is_loop =
true;
733 bool dont_prune = (options.mode & SkeletonOptions::Prune) == 0;
734 bool preserve_topology = (options.mode & SkeletonOptions::PreserveTopology) != 0 ||
735 options.mode == SkeletonOptions::Prune;
736 bool relative_pruning = (options.mode & SkeletonOptions::Relative) != 0;
737 WeightType threshold = (options.mode == SkeletonOptions::PruneTopology ||
738 options.mode == SkeletonOptions::Prune)
741 ? options.pruning_threshold*skeleton[center].salience
742 : options.pruning_threshold;
744 int branch_count = 0;
745 double average_length = 0;
746 for(
int k=0; k < (
int)raw_skeleton.size(); ++k)
748 Node p1 = raw_skeleton[k];
749 SNode & n1 = skeleton[p1];
751 if(n1.principal_child == lemon::INVALID &&
752 n1.salience >= threshold &&
756 average_length += n1.length;
757 total_length += n1.length;
760 dest[p1] = n1.salience;
761 else if(preserve_topology)
763 if(!n1.is_loop && n1.salience < threshold)
766 else if(p1 != center && n1.salience < threshold)
770 average_length /= branch_count;
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());
786 if(options.mode == SkeletonOptions::Prune)
787 detail::skeletonThinning(squared_distance, dest,
false);
790 class SkeletonFeatures
793 double diameter, total_length, average_length, euclidean_diameter;
794 UInt32 branch_count, hole_count;
795 Shape2 center, terminal1, terminal2;
801 , euclidean_diameter(0)
948 template <
class T1,
class S1,
952 MultiArrayView<2, T2, S2> dest,
953 SkeletonOptions
const & options = SkeletonOptions())
955 skeletonizeImageImpl(labels, dest, (ArrayVector<SkeletonFeatures>*)0, options);
958 template <
class T,
class S>
960 extractSkeletonFeatures(MultiArrayView<2, T, S>
const & labels,
961 ArrayVector<SkeletonFeatures> & features,
962 SkeletonOptions
const & options = SkeletonOptions())
964 MultiArray<2, float> skeleton(labels.shape());
965 skeletonizeImageImpl(labels, skeleton, &features, options);
972 #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
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
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:1775
SkeletonOptions & returnSalience()
Don'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
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
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
Option object for skeletonizeImage()
Definition: skeleton.hxx:256
SkeletonOptions & returnLength()
Don't prune and return the length of each skeleton segment.
Definition: skeleton.hxx:303
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
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
SkeletonOptions()
construct with default settings
Definition: skeleton.hxx:280
bool operator>(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
greater
Definition: fixedpoint.hxx:530
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'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