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,
398 ArrayLike * features,
401 static const unsigned int N = 2;
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.");
421 using namespace multi_math;
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);
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;
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();
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;
571 pathFinder.
reRun(weights, center, lemon::INVALID, maxWeight);
573 bool compute_salience = (options.mode & SkeletonOptions::Salience) != 0;
576 for(
int k=raw_skeleton.size()-1; k >= 0; --k)
578 Node p1 = raw_skeleton[k],
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)
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);
684 poly.push_back_unsafe(p1);
685 poly.push_back_unsafe(p2);
690 poly.push_back_unsafe(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,
960 template <
class T,
class S>
967 skeletonizeImageImpl(labels, skeleton, &features, options);
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
const Node & target() const
get the target node
Definition: graph_algorithms.hxx:427
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
double arcLengthQuantile(double quantile) const
Definition: polygon.hxx:379
const difference_type & shape() const
Definition: multi_array.hxx:1596
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
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
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
Definition: polygon.hxx:78
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
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
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
const PredecessorsMap & predecessors() const
get the predecessors node map (after a call of run)
Definition: graph_algorithms.hxx:442
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
size_type size() const
Definition: array_vector.hxx:358
SkeletonOptions()
construct with default settings
Definition: skeleton.hxx:280
bool operator>(FixedPoint< IntBits1, FracBits1 > l, FixedPoint< IntBits2, FracBits2 > r)
greater
Definition: fixedpoint.hxx:530
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'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
const DiscoveryOrder & discoveryOrder() const
get an array with all visited nodes, sorted by distance from source
Definition: graph_algorithms.hxx:437
SkeletonOptions & dontPrune()
return the un-pruned skeletong
Definition: skeleton.hxx:287