Keywords

1 Introduction

The skeletal representation of the planar object is essential for many problems of computer vision and pattern recognition, image processing, computer graphics and visualization [1]. Skeletons are widely used for shape matching [2, 3], optical character recognition [4] and image retrieval [2, 5]. In the area of biomedical image processing, skeletonization methods are extensively applied to compute the central line of thin objects. For example, one can extract the skeletal graph representing the retinal blood vessels topology [6, 7]. A similar technique can be applied to segment biological neural networks [8]. One can also use skeletonization methods to segment cellular filamentous structures using microscopy images [9, 10]. Thus, fast and accurate algorithms for computing the skeleton of the geometrical objects are of high relevance.

Related Work.

Existing algorithms for computing the skeleton can be classified based on the type of processed data. For example, morphological thinning techniques are extensively used for computing the skeleton of a binary image [11,12,13]. They allow us to obtain a pixel-level representation of the thin skeleton. On the next step, such representation can be converted into a graph using the vectorization methods described [14, 15]. However, the accuracy of the skeleton is bounded by the resolution of the pixel grid. Moreover, many of these methods are not rotation-invariant [11, 12].

Other techniques are based on central line tracing. They are commonly used to segment thin-line structures on an image (e.g., axons, dendrites of neurons [16], blood vessels [17], filamentous structures [17, 18]). These methods can directly represent the skeleton as a connected graph. However, due to the iterative nature of these methods, the execution time may vary significantly.

Another class of methods allows us to compute a skeleton of an object, whose shape is represented by simple polygons. Such polygons can be either sampled directly from a vector graphics data or can be extracted from a binary image using the tracing techniques (e.g., Marching squares [19]). Methods to construct the straight skeleton using the polygon shrinking technique with O(N log N) complexity are described in papers [20, 21]. A linear complexity method for a simple polygon without holes was introduced in [22]. A more general approach for constructing the skeleton of an arbitrary object with holes employs the Voronoi diagram [23, 24], which has computational complexity O(N log N), N is a number of primitives. In comparison to the techniques above, this approach allows us to directly compute a rotation-invariant thin skeleton of an object as a graph. Moreover, one can also employ the properties of the Voronoi diagram to solve the related geometrical problems [25] (e.g., finding fast a convex hull, nearest neighbor, maximal inscribed disk). However, due to a large number of the processed simple primitives, such method can become computationally costly. Therefore, we focus on the Voronoi-based skeletonization methods and on heuristic techniques allowing us to speed up such methods by employing the shape simplification techniques.

2 Problem Statement

We assume that a planar object has G1-continuous boundaries (except for a finite number of G0-continuous points – see critical points below). The object’s boundaries are represented by a set of simple planar polygons \( { \mathcal{S}}: = \left\{ {{\mathcal{P}}^{0} ,{\mathcal{P}}^{1} , \ldots ,{\mathcal{P}}^{m} } \right\} \), where polygon \( {\mathcal{P}}^{k} \) is defined as an ordered set of its vertices \( p_{1}^{k} ,p_{2}^{k} , \ldots ,p_{{M_{k} }}^{k} \). Polygon \( {\mathcal{P}}_{0} \) corresponds to the outer contour of the object. \( R\,{:= \mathcal{P}}_{0} \backslash \mathop {\bigcup }\nolimits_{i = 1}^{m} {\mathcal{P}}_{i} \) defines the object’s domain.

Let’s denote the set of open line segments (LS’) corresponding to the polygon \( {\mathcal{P}}_{k} \) by \( {\mathcal{L}}_{k} : = {\mathcal{L}}\left( {{\mathcal{P}}_{k} } \right) = \left\{ {l_{i}^{k} \text{ := }\left( {p_{i}^{k} ,p_{i + 1}^{k} } \right) \left| { i = 1, \ldots ,M_{k} , } \right. p_{{M_{k} + 1}}^{k} = p_{1}^{k} } \right\} \) and the set of all vertices (line segment’s endpoints) by \( {\mathcal{Q}} = \bigcup\nolimits_{k = 0}^{m} {\bigcup\nolimits_{j = 1}^{{M_{k} }} {\left\{ {p_{i}^{k} } \right\}} } \). A set \( {\mathcal{L} := }\bigcup\nolimits_{k = 1}^{m} {{\mathcal{L}}_{k} } \) contains all line segments of \( {\mathcal{S}} \).

Definition 1.

The Voronoi cell [29] corresponding to an element \( u \in {\mathcal{L}} \cup {\mathcal{Q}} \) is defined as a locus of points:

$$ {\mathcal{V}\mathcal{C}}\left( u \right) = \left\{ {p \in {\mathbb{R}}^{2} \left| {dist\left( {p,u} \right) \le dist\left( {p,w} \right),w \ne u, w \in {\mathcal{L}} \cup {\mathcal{Q}}} \right.} \right\} $$
(1)

Definition 2.

The Voronoi diagram [29] of a set of line segments \( {\mathcal{L}} \) (with endpoints \( {\mathcal{Q}} \)) is defined as a set of all Voronoi cells:

$$ {\mathcal{V}\mathcal{D}}\left( {{\mathcal{L}},{\mathcal{Q}}} \right) = \mathop {\bigcup }\limits_{{u \in {\mathcal{L}} \cup {\mathcal{Q}}}} \left\{ {{\mathcal{V}\mathcal{C}}\left( u \right)} \right\} $$
(2)

Remark 1.

The most of the computational algorithms (e.g., “Divide and Conquer” [26], Fortune’s algorithm [27]) represent the boundaries between neighboring Voronoi cells in terms of the Voronoi graph [28, 29] \( G_{{\mathcal{S}}} = \left( {V_{{\mathcal{S}}} ,E_{{\mathcal{S}}} } \right) \) with a set of the Voronoi vertices \( V_{{\mathcal{S}}} \) and a set of Voronoi edges \( E_{{\mathcal{S}}} \, \subseteq \,V_{{\mathcal{S}}} \, \times \,V_{{\mathcal{S}}} \).

Definition 3.

Let’s assume that a polygon \( {\mathcal{P}} \) approximates boundary of geometrical object and a vertex \( p \) of \( {\mathcal{P}} \) corresponds to the point of the boundary, where object is G0-continouous (but not G1-continouous). Then vertex \( p \) is called critical points (vertices) of the polygon \( {\mathcal{P}} \).

Remark 2.

Vertices of polygon \( {\mathcal{P}} \) corresponding to G1-continouous part of object’s boundary, might induce redundant edges of the Voronoi diagram – the bisectors between consecutive line segments \( l_{i} \) and \( l_{i + 1} \) sharing common non-critical endpoint \( p_{i} \). In order to obtain an approximate Voronoi diagram of an object represented by \( {\mathcal{S}} \), such redundant edges corresponding to all non-critical points of \( {\mathcal{S}} \) should be removed [29].

Definition 4.

An approximate Voronoi diagram \( {\mathcal{V}\mathcal{D}}_{a} \left( {\mathcal{S}} \right) \) [29] for a planar object represented by a set of polygons \( {\mathcal{S}} \) is obtained as a subgraph \( G_{{\mathcal{S}}}^{a} \) of the Voronoi graph \( G_{{\mathcal{S}}} \) by removing the edges of \( G_{{\mathcal{S}}} \) corresponding to the bisectors between two consecutive line segments \( l_{i} \) and \( l_{i + 1} \) sharing a common non-critical vertex \( p_{i} \).

Definition 5.

The Voronoi skeleton [23] of a planar object represented by \( {\mathcal{S}} \) is a subset of the approximate Voronoi diagram \( {\mathcal{V}\mathcal{D}}_{a} \left( {\mathcal{S}} \right) \) located inside object’s region \( R \).

Remark 3.

Thus, the Voronoi skeleton of \( {\mathcal{S}} \) is obtained by removing (or trimming) the edges of \( G_{{\mathcal{S}}}^{a} \), which do not locate in \( R \).

Problem Statement:

Given a set of polygons \( {\mathcal{S}} \), which represent a planar object, construct the Voronoi skeleton of \( {\mathcal{S}} \).

3 Algorithm

In this section, we describe the algorithm for computing the Voronoi skeleton. In Subsect. 3.2 we show an algorithm for transforming the Voronoi graph \( G_{{\mathcal{S}}} \) into the final Voronoi skeleton. The complexity analysis of the algorithm is shown in Subsect. 3.3.

3.1 Algorithm Description

Input:

\( {\mathcal{S}}: = \left\{ {{\mathcal{P}}^{1} ,{\mathcal{P}}^{2} , \ldots ,{\mathcal{P}}^{m} } \right\} \) – the set of polygons, each vertex \( p_{i}^{k} \) of the polygon \( {\mathcal{P}}^{k} \) has a binary attribute \( {\text{is}}\,{\text{Critical}}\left[ {p_{i}^{k} } \right] \in \left\{ {{\text{True}},{\text{False}}} \right\} \). Polygon \( {\mathcal{P}}^{k} \) is oriented such that its interior of the object is to the right for any its line segment (LS) \( l \in {\mathcal{P}}^{k} \).

Algorithm:

  1. 1.

    Compute Voronoi diagram of line segments \( {\mathcal{L}} \) (with endpoints \( {\mathcal{Q}} \)) \( \Rightarrow \) Obtain Voronoi graph \( G_{{\mathcal{S}}} = \left( {V_{{\mathcal{S}}} ,E_{{\mathcal{S}}} } \right) \) represented as doubly-connected edge list (DCEL) [28];

  2. 2.

    Using the breadth-first search (BFS) algorithm to traverse the Voronoi graph \( G_{{\mathcal{S}}} \) and label its edges and vertices (see Subsect. 3.2);

  3. 3.

    Remove the edges \( G_{{\mathcal{S}}} \) with labels “R” or “O” and vertices with labels “B” or “O”;

  4. 4.

    Remove isolated vertices of \( G_{{\mathcal{S}}} \) if any exist;

3.2 Labeling Voronoi Graph

We traverse the edges and vertices of the Voronoi graph \( G_{{\mathcal{S}}} \) and label them according to their role in a resulting graph of the Voronoi skeleton.

figure a

Definition 6.

Voronoi vertex (cf., Fig. 1) is called (label abbreviation is in parenthesis):

Fig. 1.
figure 1

Examples of the labeled Voronoi vertices and edges.

Full size image

  • Inner (“I”), if the vertex is located inside the object’s polygon;

  • Outer (“O”), if the vertex is located outside the object’s polygon;

  • Critical (“C”), if it coincides with one of the critical vertices of the object’s polygon;

  • Boundary (“B”), if it coincides with one of the non-critical vertices of the polygon;

Definition 7.

Voronoi edge (cf., Fig. 1) \( e \) is called (label abbreviation is in parenthesis):

  • Inner (“I”), if it locates in \( R \) and doesn’t touch (intersect) any polygon from \( {\mathcal{S}} \);

    \( \Rightarrow \) both vertices of \( e \) are labeled as Inner;

  • Critical (“C”), if it locates in \( R \) and adjacent to a critical vertex;

  • Outer (“O”), if the edge locates outside \( R \) \( \Rightarrow \) both vertices of \( e \) are labeled as Outer;

  • Redundant (“R”), if it locates in \( R \) and touches polygon’s non-critical vertex;

The pseudocode illustrating the labeling procedure and the related functions is shown in Listing 1. Firstly, we initialize the queue Q of the breadth-first search (BFS) algorithm (see function InitQueue(Q) ) by all infinite edges of the Voronoi graph \( G_{{\mathcal{S}}} \). A common data structure Label[•] is used to store labels of Voronoi edges and vertices according to the definitions 6 and 7. Then, starting from infinite edges we label all remaining edges and vertices of \( G_{{\mathcal{S}}} \) in function TraverseBFS(Q) . At each iteration of BFS algorithm we label current edge and the adjacent non-labeled vertex. In Listing 1 (Condition) ? Value1 : Value2 denotes to the ternary conditional operator.

3.3 Complexity Analysis

Lemma 1.

The complexity of Step 1 of the skeletonizing algorithm is O(N log N), where N is a number of points in a polygon.

Proof.

At the Step 1 we construct the Voronoi diagram for polygon’s line segments using Fortune’s algorithm. According to [27] the complexity of this step is O(M log M), M - number of line segments. Since \( N\sim M \), Step 1 has complexity O(N log N). ■

Lemma 2.

The complexity of Step 2 of the skeletonizing algorithm is O(N), where N is a number of the points in an input polygon.

Proof.

Step 2 is about labeling the edges and vertices of the Voronoi graph using BFS traverse algorithm. Note that the Voronoi graph is a planar connected graph. Therefore, Euler’s formula \( \left| V \right| - \left| E \right| + f = 2 \) take place, where \( \left| V \right| \), \( \left| E \right| \), \( f \) is a number of vertices, edges and faces of a graph. If \( \left| V \right| = N \), then the number of edges \( \left| E \right| = O\left( N \right) \). The BFS algorithm traverses all edges of the Voronoi graph. Since all operations within one BFS iteration can be performed in O(1), the complexity of BFS routine is \( {\text{O}}(\left| E \right| + \left| V \right|) = {\text{ O}}\left( {\text{N}} \right) \). Thus, the complexity of Step 2 is O(N). ■

Lemma 3.

The complexity of Steps 3–4 of the skeletonizing algorithm is O(N), where N is a number of the points in an input polygon.

Proof.

One edge can be removed from DCEL in O(1) by reassigning the pointers [25, 28]. According to Lemma 2, the number of edges \( \left| E \right| = O\left( N \right) \). Therefore, the complexity of Step 3 is \( O\left( N \right) \). A single isolated vertex can be removed from DCEL in O(1). Therefore, the complexity of Step 4 is \( O\left( N \right) \). ■

Theorem 1.

The complexity of the skeletonizing algorithm is O(N log N), where N is a number of the points in an input polygon.

Proof.

According to analysis of the complexities of each algorithm’s step provided in Lemmas 1–3, the total complexity of skeletonizing algorithm is O(N log N). ■

4 Optimization and Heuristics

We introduce an optimization heuristic allowing us to compute fast the Voronoi skeleton by reducing the number of vertices of input polygons. The main idea behind the optimization procedure is illustrated by the following lemma.

Lemma 4.

Let \( {\mathcal{P}} = \left\{ {p_{1} ,p_{2} , \ldots ,p_{N} } \right\} \) be a polygon and \( l_{i} \) denotes the line segment between points \( p_{i} \) and \( p_{i + 1} \) of a polygon \( {\mathcal{P}} \), \( i = 1, \ldots ,N \), (\( p_{N + 1} = p_{0} ) \). The polygon \( {\mathcal{P}} ' \) is obtained by subdividing line segments \( l_{i} , i = 1, \ldots N \) of a polygon \( {\mathcal{P}} \) such that line segment \( l_{i} \) is replaced by a polyline \( p_{i,1} ,p_{i,2} \ldots ,p_{{i,R_{i} }} \) of points on \( l_{i} \), \( i = 1,2, \ldots ,N \) (\( p_{i,1} = p_{i} \), \( p_{{i,R_{i} }} = p_{i + 1} \)). Then the Voronoi skeletons \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \) and \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) constructed using the skeletonizing algorithm above are equal (in terms of the Hausdorff distance between the corresponding Voronoi graphs).

Proof.

The Voronoi diagram of line segments of \( {\mathcal{P}} \) and \( {\mathcal{P}} ' \) consist of the bisectors of the following types: a bisector between two line segment’s interiors, a bisector between a line segment’s interior and an endpoint, bisector between two endpoints. Let’s consider these cases separately:

Case 1 (see Fig. 2a). The bisector between two line segment’s interiors \( l_{1} \) and \( l_{2} \) is a line segment \( l' \) [27, 28]. Let’s suppose that in \( {\mathcal{P}} ' \) line segment \( l_{2} \) remains the same and \( l_{1} \) is subdivided into two parts \( l_{1,1} \) and \( l_{1,2} \) connected by a shared endpoint \( q \). Then, the Voronoi cell corresponding to \( l_{1} \) in \( {\mathcal{V}\mathcal{D}}\left( {\mathcal{P}} \right) \) will be split into two Voronoi cells (corresponding \( l_{1,1} \) and \( l_{1,2} \)) of \( {\mathcal{V}\mathcal{D}}\left( {{\mathcal{P}} '} \right) \) by the Voronoi edge \( e \) such that \( e \) is a bisector between \( l_{1,1} \) and \( l_{1,2} \) which passes through \( q \) and is perpendicular to \( l_{1} \) (and therefore, \( l_{1,1} \) and \( l_{1,2} \)). Thus, the Voronoi edge \( e \) will divide bisector line segment \( l' \) in \( {\mathcal{V}\mathcal{D}}\left( {\mathcal{P}} \right) \) into two parts \( l'_{1} \) and \( l'_{2} \) in \( {\mathcal{V}\mathcal{D}}\left( {{\mathcal{P}} '} \right) \) such that \( l'_{1} \) is a Voronoi edge of the Voronoi cell of \( l_{1,1} \) and \( l'_{2} \) is a Voronoi edge of the Voronoi cell of \( l_{1,2} \). Note that \( l'_{1} \), \( l'_{2} \) and edge \( e \) are connected together by a newly introduced Voronoi vertex \( v' \). The remaining part of the Voronoi diagrams for \( {\mathcal{P}} ' \) and \( {\mathcal{P}} \) stays the same. The BFS labeling procedure (see Step 2 of the algorithm above) for Voronoi edges and vertices of \( {\mathcal{V}\mathcal{D}}\left( {{\mathcal{P}} '} \right) \) will split the introduced in \( {\mathcal{V}\mathcal{D}}\left( {{\mathcal{P}} '} \right) \) Voronoi edge \( e \) into two parts \( e_{1} \) and \( e_{2} \): one part will be labeled as “Outer” and the other part will be labeled as “Redundant”. Therefore, both parts will be removed at Step 3 of the skeletonizing algorithm and the resulting Voronoi skeleton \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) will contain the line segment edges \( l'_{1} \), \( l'_{2} \) connected by \( v' \).

Fig. 2.
figure 2

Illustration of the proof of Lemma 4

Full size image

Case 2. In case of a line segment’s interior \( l \) and an endpoint \( p \), two possible scenarios take place. First scenario is when \( p \) is an endpoint of \( l \). In this case Voronoi diagram contains an edge \( e' \) coming through \( p \) and perpendicular \( l \). The edge \( e' \) can be either removed or not by BFS procedure depending on the type of \( p \). Subdividing \( l \) into two parts \( l_{a} \) and \( l_{b} \) which share an endpoint \( q \) will introduce a new edge \( e \) parallel to \( e' \), which will be classifies as “Redundant” and removed from the final skeleton. The second scenario (see Fig. 2b) is when \( p \) is not an endpoint of \( l \). Then the bisector between \( p \) and \( l \) is a parabolic arc \( l_{p} \), which is subdivided into two parts \( l_{p,1} \), \( l_{p,2} \) if we split \( l \) into \( l_{a} \) and \( l_{b} \). The analysis in this case is the similar to the Case 1 except that now \( l'_{1} \) and \( l'_{2} \) are parabolic arcs \( l_{p,1} \) and \( l_{p,2} \), respectively.

Case 3. The bisector between two different endpoints of \( {\mathcal{V}\mathcal{D}}\left( {{\mathcal{P}} '} \right) \) or \( {\mathcal{V}\mathcal{D}}\left( {\mathcal{P}} \right) \) is an infinite edge (ray), which is classified at Step 2 of the algorithm above as “Outer” and, therefore, removed from both \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \) and \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) at Step 3.

The case of single subdivision \( (L = 1) \) of polygon’s line segment for different possible bisectors of the Voronoi diagram is covered above. The general case for several subdivisions \( L \) can be proved by induction on L.

Let’s assume that for \( L = n \) subdivisions of \( {\mathcal{P}} \) holds that \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \) and \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) are equal. The polygon \( {\mathcal{P}^{\prime\prime}} \) is obtained from \( {\mathcal{P}} ' \) by subdividing an arbitrary line segment of \( {\mathcal{P}} ' \) into two line segments. Therefore, we can apply one of the proved cases for a single subdivision above and obtain that Voronoi skeletons \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \) and \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}^{\prime\prime}}} \right) \) are equal. Thus, by induction \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \) and \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) are equal for any \( L > 0 \). ■

Remark.

It follows from Lemma 4 that the Voronoi skeleton \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) for a subdivided polygon \( {\mathcal{P}} ' \) is the same (w.r.t. Hausdorff distance) as the Voronoi skeleton \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \) for the original polygon \( {\mathcal{P}} \) (see Fig. 3). However, in comparison to \( {\mathcal{V}\mathcal{S}}\left( {\mathcal{P}} \right) \), \( {\mathcal{V}\mathcal{S}}\left( {{\mathcal{P}} '} \right) \) is represented with a larger number of Voronoi edges and vertices. Therefore, the concept of the Voronoi skeleton with a minimal number of vertices/edges take place. Applying Lemma 4 in the reverse direction allows us to reduce the number of vertices and edges of the Voronoi skeleton. This in turn reduces the execution time of skeletonization algorithm and compresses the resulting graph representation of a skeleton.

Fig. 3.
figure 3

The Voronoi skeletons (red) for polygon P (blue) and its subdivided version P’ (blue) and respective Voronoi diagrams (gray). (Color figure online)

Full size image

Therefore, our aim is to design a heuristic based on simplification operation (reverse to subdivision) and obtain polygon \( {\mathcal{P}} \) from \( {\mathcal{P}} ' \). According to Lemma 4 simplification procedure (algorithm) should meet the following requirement:

Simplification Requirement (SR):

The polygon simplification heuristic removes the points corresponding to colinear connected line segments of the polygon representing such line segments by a single line segment.

Thus, we introduce Step 0 of the skeletonizing algorithm: simplify each polygon of a set \( {\mathcal{S}} \) by reducing the points associated with colinear connected line segments (SR). This operation can be performed using one of the existing polygon simplification algorithms satisfying the simplification requirement (SR).

Analysis of Simplification Algorithms.

We have analyzed the most commonly used algorithms for polygon (polyline) simplification and summarized the results in Table 1.

Table 1. The overview of polygon (polyline) simplification algorithms.
Full size table

However, certain simplification strategies do not agree with the simplification requirement (SR). For example, a naive Nth point [38] method merely removes every Nth point from a polygon ignoring its geometry. Circle simplification [38] method groups together points forming spatial clusters based on the distance threshold. Then, a single representative point replaces each such cluster. Li-Openshaw [37] and Rapso [36] algorithms simplify polyline based on spatial pixel (or hexagon-based) grid. These algorithms instead solve the problem of polyline digitization (useful for solving the problem of optimal map rescaling). Therefore, we consider only the algorithms fulfilling SR (see Table 2). Most of the analyzed algorithms have complexity O(N) except DP [30] and VW [31] algorithms with O(N log N) complexity. In order to select the algorithm, which shows the best performance improvement, has the smallest computational overhead and influences the resulting skeleton the least, we investigated these algorithms empirically as described in the evaluation section.

Table 2. Suitable polygon simplification algorithms, their parameter and heuristics.
Full size table

5 Evaluation

We evaluate the performance of the skeletonization algorithm in terms of the execution time and measure the influence of the introduced heuristics onto the accuracy, execution time of the overall algorithm. We also estimate the computational overheads related to the line simplification algorithms.

Dataset.

In order to evaluate the performance of the skeletonization algorithm and individual optimization heuristics, we used polygons obtained from MPEG 7 CE-Shape-1 dataset. These polygons were extracted from binary images using the Marching Squares algorithm [19]. In total the dataset consists of 1282 polygons (see Fig. 4).

Fig. 4.
figure 4

Distribution of polygon’s sizes

Full size image

Measures.

We have measured the following quantities:

  1. 1.

    Execution time (ms) of each simplification algorithm, skeletonizing algorithm with (without) the mentioned heuristics and overall execution time. The experiments were carried on Intel Core i7, 2.2 GHz, 16 Gb RAM.

  2. 2.

    Hausdorff distances \( d_{H} \) (errors) [39] between the simplified and original polygons and also between the ground truth skeleton and one obtained using the skeletonization with heuristics;

  3. 3.

    Simplification rate (%) of the polygon is computed as follows:

$$ SR\left( {P,P^{\prime}} \right) = \frac{{\left| P \right| - \left| {P'} \right|}}{\left| P \right|} \cdot 100\% $$
(3)

where \( P \) and \( P^{\prime} \) are original and simplified polygons, respectively. \( \left| P \right| \) is the number of vertices of \( P \) (large values of \( SR\left( {P,P^{\prime}} \right) \) correspond to small \( \left| {P'} \right| \) w.r.t. \( \left| P \right| \)).

Parameters.

The parameters of the simplification algorithms (see Table 2) were chosen using the line search method such that the maximum simplification rate is achieved for a given threshold value of Hausdorff distance \( d_{H} \) between the simplified polygon \( P^{\prime} \) and an original polygon \( P \). This allows us to compare different simplification algorithm with respect to the maximum tolerable error. The established parameters of the simplification algorithms for the respective values of \( d_{H} \) are shown in Table 3.

Table 3. Parameters of the simplification algorithms.
Full size table

For the algorithms with two parameters we applied additional heuristics to choose the value of the second parameter (see Table 2). These heuristics were devised to achieve the maximum simplification rate for a given Hausdorff error threshold \( d_{H} \). It was established that for LA and PD algorithms the optimal value of \( \theta \) is \( 0.25 \) (for \( \theta > 0.25 \) the simplification rate does not increase, but the execution time of these simplification algorithms rises).

Evaluation Results.

We have measured the execution time of each suitable simplification algorithm for fixed values of Hausdorff error thresholds \( d_{H} \) (see Fig. 5a). These measurements show the computational overheads related to the optimization step of the skeleton algorithm. In order to compare the quality of the simplification algorithms, we measured the respective simplification rates for given values of \( d_{H} \).

Fig. 5.
figure 5

Execution time (ms), simplification rates (%) of optimization heuristics

Full size image

Figure 5b shows that the algorithms of LA and ZS have the most substantial extent of polygon simplification (compression) for a given \( d_{H} \) having nearly identical dependency curves. PD, VW, and PD algorithms achieve slightly smaller simplification rates showing almost undistinguishable behavior for most of the cases. However, VW algorithm overperforms other algorithms for small values of \( d_{H} \, < \,0.002 \). OP and RW algorithms have the lowest simplification rates with nearly identical dependency curves.

In Fig. 5 one notices that despite being the fastest, algorithms of OP and RW have the smallest simplification rate and, therefore, might not guarantee the fastest execution of the skeletonization algorithm. Therefore, we measured the total execution time of the skeletonization algorithm depending on the value of \( d_{H} \) taking into account the overhead time of the simplification heuristics (see Fig. 6a).

Fig. 6.
figure 6

Total average execution time (ms) and skeletonization errors;

Full size image

Based on Fig. 6a we can choose the fastest optimization heuristics. However, different values of \( d_{H} \) threshold might affect the accuracy of the final skeleton. Therefore, we investigated the influence of \( d_{H} \) on the result of the skeletonization algorithm. We calculated the skeletonization error as Hausdorff distance between the ground truth skeleton and the result of optimized skeletonization algorithm (see Fig. 6b).

6 Discussion

Figure 6a shows that DP- and VW-based heuristics reduce the computational time to the greatest extent. Only these two heuristics overperform the optimization-free approach (NO) for small values of \( d_{H} \le 0.001 \). The optimization based on OP and RW algorithms shows the smallest skeletonization error among the other approaches (see Fig. 6b). However, for \( d_{H} \, < \,0.002 \) these algorithms have outsized computational overheads eliminating the whole effect of the optimization. Therefore, it is reasonable to use OP and RW algorithm only for \( d_{H} \, > \,0.002 \). Note that the variance of skeletonization error for different heuristics decreases as \( d_{H} \to 0 \) (see Fig. 6b).

We computed 2-sample t-test to validate the hypothesis that DP- and VW-based optimizations produce different average skeletonization errors. The test showed that the errors produced by DP and VW optimizations are undistinguishable (p-value \( \approx 0.24 \)).

Another hypothesis testing was performed to distinguish the execution time between DP and VW heuristics. It showed that for the most of the cases (except \( d_{H} = 0.001 \)) DP algorithm overperforms VW (p-value \( < 0.001 \)).

Speed-Accuracy Trade-Off.

Figure 6 shows that none of the tested algorithms minimizes the accuracy and execution time of the skeletonizing method at the same time. Therefore, the choice of the heuristics is a trade-off between accuracy and the execution time. Based on the performed computational experiments the following conclusions are drawn:

  1. 1.

    If accuracy of the resulting skeleton is critical, then for \( d_{H} > 0.002 \) the optimization can be performed using OP or RW algorithms. However, for \( d_{H} < 0.002 \) the only reasonable optimization is using the DP or VW algorithms;

  2. 2.

    If execution time of the algorithm is more critical than the accuracy, then optimization can be performed using DP or VW algorithms, which according to the provided experiments give 1.7 times less accurate result then RW and OP heuristics;

Pruning Effect of Polygon Simplification.

It was experimentally discovered, that the introduced optimization heuristics influences the skeleton in a similar way as pruning methods [40]. Figure 7 shows that for large values of \( d_{H} \) (see bottom row) simplification heuristics tends to regularize shape of the object in a way that the branches of the skeleton corresponding to small shape perturbation disappear (cf., Fig. 7, top row). Therefore, such optimization allows us not only to speed-up the execution of the skeletonization, but also to achieve a pruning effect and remove the noisy branches of the skeleton.

Fig. 7.
figure 7

Examples of the optimized Voronoi skeletons for shapes from MPEG 7 CE-Shape-1 dataset. Optimization heuristics is DP. For \( {\text{d}}_{\text{H}} = 0.001 \) (top row of images) skeletons contain redundant branches in comparison to \( {\text{d}}_{\text{H}} = 1.0 \) (the bottom row of images).

Full size image

7 Conclusion

We proposed optimization heuristics for computing the Voronoi skeleton of the polygonal data. This topic is of relevance due to its direct relation to efficient processing of the vectorized geometrical representations (e.g., for image processing, computer vision, computer graphics). We illustrated in detail the main steps of the Voronoi-based skeletonization algorithm and determined that its complexity is O(N log N), where N is the number of vertices in a polygon. We also established an optimization criterion (requirement) and proposed theoretically justified optimization heuristics based on the polygon simplification algorithms. In order to evaluate the efficiency of the proposed heuristic, a series of computational experiments were conducted using the polygons from MPEG 7 CE-Shape-1 dataset. Seven state-of-the-art simplification algorithms were evaluated to determine the most suitable optimization heuristic fulfilling the established criterion. We measured the execution time of the skeletonization algorithm with and without the heuristic optimizations and determined the computational overheads related to such heuristics. We also determined the accuracy of the skeleton produced by the optimized algorithm based on the proposed heuristics. As a result, we established the criteria, which allow us to choose the optimal heuristics depending on the system’s requirement. For example, DP- and VW-based heuristics allow us to speed up the skeleton computation at least by 30%. It was discovered experimentally, that the optimization heuristics have a pruning effect onto the resulting skeleton.