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 weights[*arc] = infiniteWeight;
680 pathFinder.
reRun(weights, p1, p2);
683 poly.push_back_unsafe(p1);
684 poly.push_back_unsafe(p2);
689 poly.push_back_unsafe(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,
958 template <
class T,
class S>
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
Define a grid graph in arbitrary dimensions.
Definition: multi_fwd.hxx:217
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'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:2474
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: 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:206
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:1648
Option object for skeletonizeImage()
Definition: skeleton.hxx:256
const Node & target() const
get the target node
Definition: graph_algorithms.hxx:427
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
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:704
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'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