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"
56 Node parent, principal_child;
57 double length, salience;
62 : parent(lemon::INVALID)
63 , principal_child(lemon::INVALID)
70 SkeletonNode(Node
const & s)
72 , principal_child(lemon::INVALID)
83 typedef SkeletonNode<Node> SNode;
84 typedef std::map<Node, SNode> Skeleton;
86 Node anchor, lower, upper;
90 : anchor(lemon::INVALID)
95 void addNode(Node
const & n)
97 skeleton[n] = SNode(n);
99 lower = min(lower, n);
100 upper = max(upper, n);
104template <
class Graph,
class Node,
class NodeMap>
106neighborhoodConfiguration(Graph
const & g, Node
const & n, NodeMap
const & labels)
108 typedef typename Graph::OutArcIt ArcIt;
109 typedef typename NodeMap::value_type LabelType;
111 LabelType label = labels[n];
113 for(ArcIt arc(g, n); arc != lemon::INVALID; ++arc)
115 v = (v << 1) | (labels[g.target(*arc)] == label ? 1 : 0);
120template <
class Node,
class Cost>
121struct SkeletonSimplePoint
126 SkeletonSimplePoint(Node
const & p, Cost c)
130 bool operator<(SkeletonSimplePoint
const & o)
const
132 return cost < o.cost;
135 bool operator>(SkeletonSimplePoint
const & o)
const
137 return cost > o.cost;
141template <
class CostMap,
class LabelMap>
143skeletonThinning(CostMap
const & cost, LabelMap & labels,
144 bool preserve_endpoints=
true)
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;
152 typedef SkeletonSimplePoint<Node, double> SP;
155 std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue;
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,
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
185 bool * isSimplePoint = preserve_endpoints
186 ? isSimplePreserveEndPoints
189 int max_degree = g.maxDegree();
190 double epsilon = 0.5/labels.size(), offset = 0;
191 for (NodeIt node(g); node != lemon::INVALID; ++node)
194 if(g.out_degree(p) == max_degree &&
196 isSimplePoint[neighborhoodConfiguration(g, p, labels)])
199 pqueue.push(SP(p, cost[p]+offset));
206 Node p = pqueue.top().point;
210 !isSimplePoint[neighborhoodConfiguration(g, p, labels)])
217 for (ArcIt arc(g, p); arc != lemon::INVALID; ++arc)
219 Node q = g.target(*arc);
220 if(g.out_degree(q) == max_degree &&
222 isSimplePoint[neighborhoodConfiguration(g, q, labels)])
224 pqueue.push(SP(q, cost[q]+offset));
231template <
class Label,
class Labels>
235 Labels
const & labels;
237 CheckForHole(Label l, Labels
const & ls)
242 template <
class Node>
243 bool operator()(Node
const & n)
const
245 return labels[n] == label;
263 PreserveTopology = 4,
266 PruneCenterLine = 32,
267 PruneLength = Length + Prune,
268 PruneLengthRelative = PruneLength + Relative,
269 PruneSalience = Salience + Prune,
270 PruneSalienceRelative = PruneSalience + Relative,
271 PruneTopology = PreserveTopology + Prune
275 double pruning_threshold;
282 : mode(SkeletonMode(PruneSalienceRelative | PreserveTopology))
283 , pruning_threshold(0.2)
298 mode = PruneCenterLine;
319 if(preserve_topology)
320 mode = SkeletonMode(mode | PreserveTopology);
321 pruning_threshold = threshold;
332 mode = PruneLengthRelative;
333 if(preserve_topology)
334 mode = SkeletonMode(mode | PreserveTopology);
335 pruning_threshold = threshold;
355 mode = PruneSalience;
356 if(preserve_topology)
357 mode = SkeletonMode(mode | PreserveTopology);
358 pruning_threshold = threshold;
369 mode = PruneSalienceRelative;
370 if(preserve_topology)
371 mode = SkeletonMode(mode | PreserveTopology);
372 pruning_threshold = threshold;
386 mode = PruneTopology;
393template <
class T1,
class S1,
397skeletonizeImageImpl(MultiArrayView<2, T1, S1>
const & labels,
398 MultiArrayView<2, T2, S2> dest,
399 ArrayLike * features,
400 SkeletonOptions
const & options)
402 static const unsigned int N = 2;
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;
414 vigra_precondition(labels.shape() == dest.shape(),
415 "skeleton(): shape mismatch between input and output.");
417 MultiArray<N, MultiArrayIndex> squared_distance;
422 using namespace multi_math;
424 MultiArray<N, Shape> vectors(labels.shape());
428 ArrayVector<Node> ends_to_be_checked;
429 Graph g(labels.shape());
430 for (EdgeIt edge(g); edge != lemon::INVALID; ++edge)
432 Node p1 = g.u(*edge),
436 maxLabel = max(maxLabel, max(l1, l2));
442 const Node v1 = vectors[p1],
446 if(max(
abs(dv)) <= 1)
456 if(squared_distance[p1] == 4)
457 ends_to_be_checked.push_back(p1);
462 if(squared_distance[p2] == 4)
463 ends_to_be_checked.push_back(p2);
468 if(l1 > 0 && max(
abs(vectors[p1] + p1 - p2)) > 1)
470 if(l2 > 0 && max(
abs(vectors[p2] + p2 - p1)) > 1)
479 for (
unsigned k=0; k<ends_to_be_checked.size(); ++k)
483 Node p1 = ends_to_be_checked[k];
486 for (ArcIt arc(g8, p1); arc != lemon::INVALID; ++arc)
488 Node p2 = g8.target(*arc);
489 if(dest[p2] == label && squared_distance[p2] < 4)
493 dest[p1+vectors[p1]/2] = label;
506 detail::skeletonThinning(squared_distance, dest);
508 if(options.mode == SkeletonOptions::PruneCenterLine)
514 features->resize((
size_t)maxLabel + 1);
515 ArrayVector<detail::SkeletonRegion<Node> > regions((
size_t)maxLabel + 1);
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)
528 regions[(size_t)label].addNode(p1);
530 for (ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
532 Node p2 = g.target(*arc);
533 if(dest[p2] == label)
534 weights[*arc] =
norm(p1-p2);
536 weights[*arc] = infiniteWeight;
540 ShortestPathDijkstra<Graph, WeightType> pathFinder(g);
542 for(std::size_t label=1; label < regions.size(); ++label)
544 Skeleton & skeleton = regions[label].skeleton;
545 if(skeleton.size() == 0)
549 Node anchor = regions[label].anchor,
550 lower = regions[label].lower,
551 upper = regions[label].upper + Shape(1);
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();
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()]);
563 if(options.mode == SkeletonOptions::PruneCenterLine)
565 for(
unsigned int k=0; k<center_line.size(); ++k)
566 dest[center_line[k]] = (T2)label;
571 Node center = center_line[
roundi(center_line.arcLengthQuantile(0.5))];
572 pathFinder.reRun(weights, center, lemon::INVALID, maxWeight);
574 bool compute_salience = (options.mode & SkeletonOptions::Salience) != 0;
575 ArrayVector<Node> raw_skeleton(pathFinder.discoveryOrder());
577 for(
int k=raw_skeleton.size()-1; k >= 0; --k)
579 Node p1 = raw_skeleton[k],
580 p2 = pathFinder.predecessors()[p1];
581 SNode & n1 = skeleton[p1];
582 SNode & n2 = skeleton[p2];
587 for (BackArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
589 Node p = g.target(*arc);
590 if(weights[*arc] == infiniteWeight)
594 if(pathFinder.predecessors()[p] == p1)
596 if(n1.principal_child == lemon::INVALID ||
597 skeleton[p].principal_child == lemon::INVALID)
599 weights[*arc] = infiniteWeight;
603 WeightType l = n1.length +
norm(p1-p2);
607 n2.principal_child = p1;
612 const double min_length = 4.0;
614 if(n1.length >= min_length)
616 n1.salience = max(n1.salience, (n1.length + 0.5) / sqrt(squared_distance[p1]));
619 if(n2.salience < n1.salience)
620 n2.salience = n1.salience;
623 else if(options.mode == SkeletonOptions::DontPrune)
624 n1.salience = dest[p1];
626 n1.salience = n1.length;
630 for(
int k=0; k < (int)raw_skeleton.size(); ++k)
632 Node p1 = raw_skeleton[k];
633 SNode & n1 = skeleton[p1];
635 SNode & n2 = skeleton[p2];
637 if(p1 == n2.principal_child)
639 n1.length = n2.length;
640 n1.salience = n2.salience;
644 n1.length +=
norm(p2-p1);
646 n1.partial_area = n2.partial_area + (p1[0]*p2[1] - p1[1]*p2[0]);
651 skeleton[center].is_loop =
true;
655 detail::CheckForHole<std::size_t, MultiArrayView<2, T1, S1> > hasNoHole(label, labels);
657 double total_length = 0.0;
658 for(
int k=raw_skeleton.size()-1; k >= 0; --k)
660 Node p1 = raw_skeleton[k];
661 SNode & n1 = skeleton[p1];
663 if(n1.principal_child == lemon::INVALID)
665 for (ArcIt arc(g, p1); arc != lemon::INVALID; ++arc)
667 Node p2 = g.target(*arc);
668 SNode * n2 = &skeleton[p2];
672 if(weights[*arc] == infiniteWeight)
675 MultiArrayIndex area2 =
abs(n1.partial_area - (p1[0]*p2[1] - p1[1]*p2[0]) - n2->partial_area);
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(
decltype(poly.size()) 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];
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);
790class 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)
948template <
class T1,
class S1,
958template <
class T,
class S>
960extractSkeletonFeatures(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);
Definition array_vector.hxx:514
TinyVector< MultiArrayIndex, N > type
Definition multi_shape.hxx:272
Base class for, and view to, vigra::MultiArray.
Definition multi_fwd.hxx:127
@ OuterBoundary
Pixels just outside of each region.
Definition multi_distance.hxx:836
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition sized_int.hxx:183
void skeletonizeImage(...)
Skeletonization of all regions in a labeled 2D image.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition fixedpoint.hxx:616
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition fixedpoint.hxx:1775
@ IndirectNeighborhood
use direct and indirect neighbors
Definition multi_fwd.hxx:188
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
bool operator>(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
greater
Definition fixedpoint.hxx:530
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
bool operator<(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
less than
Definition fixedpoint.hxx:512
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
Option object for skeletonizeImage()
Definition skeleton.hxx:258
SkeletonOptions & pruneSalience(double threshold, bool preserve_topology=true)
prune skeleton segments whose salience is below the given threshold
Definition skeleton.hxx:353
SkeletonOptions & returnSalience()
Don't prune and return the salience of each skeleton segment.
Definition skeleton.hxx:341
SkeletonOptions & pruneLength(double threshold, bool preserve_topology=true)
prune skeleton segments whose length is below the given threshold
Definition skeleton.hxx:316
SkeletonOptions & returnLength()
Don't prune and return the length of each skeleton segment.
Definition skeleton.hxx:304
SkeletonOptions & pruneCenterLine()
return only the region's center line (i.e. skeleton graph diameter)
Definition skeleton.hxx:296
SkeletonOptions()
construct with default settings
Definition skeleton.hxx:281
SkeletonOptions & pruneSalienceRelative(double threshold, bool preserve_topology=true)
prune skeleton segments whose relative salience is below the given threshold
Definition skeleton.hxx:367
SkeletonOptions & pruneTopology(bool preserve_center=true)
prune such that only the topology is preserved
Definition skeleton.hxx:383
SkeletonOptions & pruneLengthRelative(double threshold, bool preserve_topology=true)
prune skeleton segments whose relative length is below the given threshold
Definition skeleton.hxx:330
SkeletonOptions & dontPrune()
return the un-pruned skeletong
Definition skeleton.hxx:288