AN OPTIMIZED ALGORITHM FOR COMPUTING THE VORONOI SKELETON

The skeleton-based representation is widely used in such fields as computer graphics, computer vision and image processing. Therefore, efficient algorithms for computing planar skeletons are of high relevance. In this paper, we propose an optimized algorithm for computing the Voronoi skeleton of a planar object with holes, which is represented by a set of polygons. Such skeleton allows us to use efficiently the properties of the underlying Voronoi diagram data structure. It was shown that the complexity of the proposed Voronoi-based skeletonization algorithm is O(N log N), where N is the number of polygon’s vertices. We have also proposed theoretically justified optimization heuristic based on polygon/polyline simplification algorithms. In order to evaluate and prove the efficiency of the heuristic, a series of computational experiments were conducted involving the polygons obtained from MPEG 7 CE-Shape-1 dataset. We have measured the execution time of the skeletonization algorithm, computational overheads related to the introduced heuristics and also the influence of the heuristic onto accuracy of the resulting skeleton. As a result, we have established the criteria, which allow us to choose the optimal heuristics for our skeletonization algorithm depending on the system’s requirements.


INTRODUCTION
The skeletal representation of the planar object is essential for many problems of computer vision and pattern recognition, computer graphics and visualization [1]. For example, skeletons are widely used for shape matching [2,3], optical character recognition [4] and image retrieval [2,5]. In the biological image processing, skeletonization methods are extensively applied to extract the central line of thin objects. For example, based on the image data one can obtain the skeletal graph representing the retinal blood vessels topology [6,7]. Similarly, this technique is applied for segmenting the biological neurons [8] and for extracting thin subcellular structures [9,10] based on microscopy data.

Related work.
There are different skeletonization methods depending on a type of input data. Morphological thinning techniques are extensively used for computing the skeleton of a binary image [11][12][13]. They allow us to obtain a pixel-based representation of the skeleton. In order to obtain the topology of the underlying graph, graph vectorization methods can be applied [14,15].
Another type of techniques is based on central line tracing. They are commonly used to segment thin-structures depicted on an image (e.g., axons and 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 an iterative nature of these methods, the execution time may vary significantly. A different type of methods is used to compute a skeleton of an object, whose shape is represented by polygonal contours (with holes). Such contours are extracted from a binary image using tracing techniques (e.g., Marching squares [19]). Alternatively, these methods can compute a skeleton of primitives obtained from a vector graphic representation.
Methods to construct the straight skeleton based on 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].
Another approach for constructing the skeleton of an object with polygonal representation is by using the Voronoi diagram [23,24]. This class of methods allows us to obtain a graph representation of the skeleton directly and use the advantages of the Voronoi diagram data structure allowing for solutions of many related problems [25] (e.g., finding a convex hull, nearest neighbor, maximal inscribed disk). However, the main drawback of such methods is high computational costs due to processing the vast number of simple primitives (points, line segments).
Therefore, in this paper, we propose an approach to decrease the computational costs related to the Voronoi-based skeletonization algorithms by introducing the novel preprocessing step based on shape simplification heuristics. Our primary focus is on investigation the suitable optimization heuristics, empirical validation of such heuristics and determining their properties.
The paper is organized as follows: The first section gives an overview of the problem and related state-of-the-art literature. In the second section, we outline the theoretical background and formulate the skeletonization problem. The third section is devoted to the description of the skeletonization algorithm, its main steps, and vertex/edge labeling procedure. In Section 4 we introduce theoretical background related to the optimization heuristics, we formulate the shape simplification requirement and analyze suitable shape simplification algorithms. In Section 5 we show the main results of our empirical study of the optimized skeletonization algorithm; we illustrate the computational overheads related to optimization heuristics and demonstrate the algorithm's performance improvements. Section 6 is devoted to the analysis of the obtained results. It establishes the criteria allowing the selection of the optimal heuristics for Voronoi-based skeletonization algorithm depending on the system's requirements. Finally, the paper concludes with a summary of the main results.

PROBLEM STATEMENT
We assume that the planar object has continuously-differentiable boundaries (called G1continuous) except for a finite number of critical points, where contour is continuous, but not differentiable (G0-continuous points the set of all line segments representing the planar object. We also assume that the polygon 0 represents the outer contour of the object and polygons 1 , … correspond to the inner holes. Thus, ≔ 0 \ ⋃ =1 defines the object's region (domain). Definition 1. The Voronoi cell corresponding to an element ∈ ℒ ∪ is a locus of points: The Voronoi diagram of a set of line segments ℒ (with endpoints ) is defined as a set of all Voronoi cells: Remark 1. Note that the Voronoi diagram of line segments ℒ (with endpoints ) is defined by the boundaries between the neighboring Voronoi cells. The most of the computational algorithms (e.g., "Divide and Conquer" [26], Fortune's algorithm [27]) represent such boundaries in terms of the Voronoi graph [28][29] = ( , ) with a set of the Voronoi vertices and a set of Voronoi edges ⊆ × . Definition 3. A finite set of polygon's points, where the object is G0-continouous (but not G1continouous) are called critical points (vertices) of the polygon .
Remark 2. Note that the vertices of the polygon , which correspond to G1-continouous part of the object's boundary, induce redundant edges of the Voronoi diagramthe bisectors between two consecutive line segments and +1 , which share a common non-critical vertex . In order to obtain an approximate Voronoi diagram of an object represented by , such redundant edges corresponding to all non-critical points of should be removed. Definition 4. The approximate Voronoi diagram ( ) for a planar object represented by a set of polygons is obtained as a subgraph of the Voronoi graph obtained by removing the edges of corresponding to the bisectors between two consecutive line segments and +1 , which share a common non-critical vertex . Definition 5. The Voronoi skeleton of a planar object represented by a set of polygons is a subset of the approximate Voronoi diagram ( ) positioning inside the object's region .
Remark 3. Thus, the Voronoi skeleton of is obtained by removing (or trimming) the edges of , which do not belong to Problem statement: Given a set of polygons , which represent a planar object, construct the Voronoi skeleton of .

SKELETONIZATION METHOD
In this section, we provide a description of the Voronoi skeletonization algorithm and Voronoi graph processing routine. We also analyze the algorithm's complexity (see Subsection 3.3).

SKELETONIZATION ALGORITHM
The input of the algorithm and its main steps are the following: Polygon is oriented such that the interior of the object is to the right for any line segment ∈ .

LABELING THE VORONOI GRAPH
We traverse the edges and vertices of the Voronoi graph and mark them according to their role in a resulting graph of the Voronoi skeleton. Definition 6. The Voronoi vertex is called: ▪ Inner, if it is inside the object's polygon; ▪ Outer, if it is outside the object's polygon; ▪ Critical, if it coincides with one of the critical vertices of the object's polygon; ▪ Boundary, if it coincides with one of the non-critical vertices of the polygon; Definition 7. The Voronoi edge is called: ▪ Inner, if it is located inside object's region and doesn't touch its boundaries; ▪ Critical, if it is located inside object's region and adjacent to a critical vertex; ▪ Outer, if the edge is located outside the object's polygon; ▪ Redundant, if it is located inside object's polygon and touches polygon's noncritical vertex; Therefore, the labels for Voronoi vertices are (see Fig. 2) Inner -"I" Critical -"C", Outer -"O", Boundary -"B"; The labels for Voronoi edges are Inner -"I", Critical -"C", Outer -"O", Redundant -"R" (see Fig. 2). The types of the Voronoi cells are Endpoint -"EP" and Line segment -"LS"; Note that Outer and Boundary vertices are absent in the final Voronoi skeleton and therefore should be removed from the Voronoi graph together with Outer and Redundant edges.
We label the edges and vertices of iteratively using the breadth-first search procedure (BFS) Firstly, we initialize the queue of BFS algorithm by infinite edges of a Voronoi graph. The pseudocode for the BFS initialization procedure is shown in the listing below:

COMPLEXITY ANALYSIS
The total complexity of the skeletonizing algorithm above is described in Theorem 1, which uses Lemmas 1-3 to establish the complexities of each individual step of the algorithm. 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. The Step 1 is about the construction of the Voronoi diagram for a line segments, which are extracted from the input polygons. The complexity of the Fortune's algorithm for Voronoi diagram construction algorithm according to [27]  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 | | − | | + = 2 takes place, where | |, | |, is a number of vertices, edges and faces of a graph. If | | = , then the number of edges | | = ( ).
The BFS algorithm traverses all edges of the Voronoi graph. Since all operations within one BFS iteration can be performed in 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 | | = ( ). Therefore, the complexity of Step 3 is ( ). A single isolated vertex can be removed from DCEL in O (1). Therefore, the complexity of Step 4 is ( ). ■ 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). ■

OPTIMIZATION
The purpose of this section is to introduce an optimization heuristic algorithm, which allow us to compute fast the Voronoi skeleton by reducing the number of vertices of input polygons. The main idea behind the optimized Voronoi skeleton construction algorithm is illustrated by the following lemma.  Proof. The Voronoi diagram of line segments of and ′ consists of the bisectors of the following types: a bisector between two line segment interiors, a bisector between a line segment interior and an endpoint and a bisector between two endpoints. Let us consider these cases separately: Case 1 (see Fig. 4). The bisector between two line segment interiors 1 and 2 is a line segment ′ [27,28]. Let us suppose that in ′ line segment 2 remains the same and 1 is subdivided into two parts 1,1 and 1,2 connected by a shared endpoint . Then, the Voronoi cell corresponding to 1 in ( ) will be split into two Voronoi cells (correspondingly 1,1 and 1,2 ) of ( ′) by the Voronoi edge such that is a bisector between 1,1 and 1,2 which passes through and is perpendicular to 1 (and therefore, 1,1 and 1,2 ). Thus, the Voronoi edge will divide bisector line segment ′ in ( ) into two parts ′ 1 and ′ 2 in ( ′) such that ′ 1 is a Voronoi edge of the Voronoi cell of 1,1 and ′ 2 is the Voronoi edge of the Voronoi cell of 1,2 . Note that ′ 1 , ′ 2 and edge are connected together by a newly introduced Voronoi vertex ′. The remaining part of the Voronoi diagrams for ′ and stays the same.
The BFS labeling procedure (see Step 2 of the algorithm above) for Voronoi edges and vertices of ( ′) will split the introduced in ( ′) Voronoi edge into two parts 1 and 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 ( ′) will contain the line segment edges ′ 1 , ′ 2 connected by ′.
Case 2. In case of a line segment interior and an endpoint , two possible scenarios take place. First scenario is when is an endpoint of . In this case Voronoi diagram contains an edge ′ coming through and perpendicular . The edge ′ can be either removed or not by BFS procedure depending on the type of . Subdividing into two parts and which share an endpoint will introduce a new edge parallel to ′, which will be classified as "Redundant" and removed from the final skeleton. The second scenario (see Fig. 5) is when is not an endpoint of . Then the bisector between and is a parabolic arc , which is subdivided into two parts ,1 , ,2 if we split into and . The analysis in this case is the similar to the Case 1 except that now ′ 1 and ′ 2 are parabolic arcs ,1 and ,2 , respectively.

Figure 6 -The Voronoi skeletons (red) for polygon (blue) and its subdivided version ′ (blue) and respective Voronoi diagrams (gray)
Case 3. The bisector between two different endpoints of ( ′) or ( ) is an infinite edge (ray), which is classified at Step 2 of the algorithm above as "Outer" and, therefore, removed from both ( ) and ( ′) at Step 3. The case of single subdivision ( = 1) of the polygon line segment for different possible bisectors of the Voronoi diagram is covered above. The general case for several subdivisions can be proved by induction on L as described below.
Let us assume that for = subdivisions of hold that ( ) and ( ′) are equal. The polygon ′′ is obtained from ′ by subdividing an arbitrary line segment of ′ into two line segments. Therefore, we can apply one of the proved cases for a single subdivision above and obtain that Voronoi skeletons ( ) and ( ′′) are equal. Thus, by induction ( ) and ( ′) are equal for any > 0. ■

Remark. It follows from Lemma 4 that the Voronoi skeleton
( ′) for a subdivided polygon ′ is the same (w.r.t. Hausdorff distance) as the Voronoi skeleton ( ) for the original polygon (see Fig. 6). However, in comparison to ( ), ( ′) is represented with a larger number of Voronoi edges and vertices. Therefore, the concept of the Voronoi skeleton with a minimal number of vertices and edges take place. By applying Lemma 4 in the reverse direction, one aims to reduce the number of vertices and edges in the Voronoi skeleton. This in turn allows us to reduce the execution time of Voronoi skeletonization algorithm and also to compress the resulting graph representation of a skeleton preserving its geometrical properties. Therefore, the operation reverse to subdivisionsimplification, should be applied to polygon ′ in order to obtain . According to Lemma 4 simplification procedure (algorithm) should meet the following requirement: Simplification requirement: The polygon simplification heuristic should reduce the points corresponding to colinear consecutive line segments of the polygon. The polylines formed by such points should be replaced by a single line segment.
Thus, we introduce the Step 0 in the skeletonizing algorithm: simplify each polygon of a set by reducing the points associated with colinear consecutive line segments of the corresponding polygon. This operation can be performed using one of the existing polygon simplification algorithms, which satisfies the simplification requirement above.
Analysis of simplification algorithms. We have analyzed the most commonly used algorithms for polygon (polyline) simplification and summarized the results in Table 1.
The aim of the mentioned in Table 1 simplification algorithms is to reduce the number of points representing the polygon (polyline). However, certain simplification strategies do not agree with the simplification requirement derived from Lemma 4. For example, a naive Nth point simplification [38] method merely removes each Nth point from a polygon ignoring its geometry. Another algorithm -Circle simplification [38], aims to group together points forming spatial clusters based on the distance threshold and replace these clusters by a single representative. Li-Openshaw [37] and Rapso [36] algorithms simplify polyline based on spatial pixel (or hexagon-based) grid. The latter two algorithms rather solve the problem of polyline digitization (useful, for example, for solving the problem of optimal map rescaling). Therefore, we considered only the algorithms fulfilling the simplification requirement above (see Table 2). Note that most algorithms in Table 1 have linear complexity (except Ramer-Douglas-Peucker [30] and Visvalingam-Whyatt [31], which have O (N log N) complexity). The open issue is to choose the simplification algorithm, which would allow us to achieve the best performance improvement showing the minimum influence on the resulting skeleton. This issue is empirically investigated in the following evaluation section.

ALGORITHM EVALUATION
In this section we evaluate the performance of skeletonization algorithm in terms of execution time and measure the influence of the heuristic optimization step onto the accuracy and execution time of the overall algorithm. We also evaluate the computational overheads related to suitable line simplification algorithms. Dataset. In order to empirically evaluate the performance of the overall skeletonization algorithm and individual heuristic simplification algorithms we used polygons obtained from the dataset MPEG 7 CE-Shape-1. These polygons were extracted from binary images using Marching Squares algorithm [19]. In total our dataset contains 1282 polygons. The distribution of polygon sizes (number of vertices) is shown in Fig. 7.
Measures. In the performed experiments we have measured the following quantities: 1. Execution time (ms) of each simplification algorithm, skeletonizing algorithm with (without) the mentioned heuristics and the total execution time. The experiments were carried out on Intel Core i7, 2.2GHz, 16Gb RAM.
2. Hausdorff distances (errors) [39] between the simplified polygon and original polygon and also between the ground truth skeleton and the skeleton obtained using the skeletonization with a simplification heuristic; 3. Simplification rate (%) of the polygon is computed as follows: where is an original polygon, ′ is a simplified polygon and | | is a number of vertices of a polygon . High simplification rate means that simplified polygon has a small number of vertices in comparison to the original one.
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 the Hausdorff error the distance between simplified polygon ′ and an original polygon . This allows us to compare different simplification algorithms with respect to the maximum allowed error. The established parameters of the simplification algorithms for the respective values of are shown in Table 3.

Figure 8 -Execution time of the simplification algorithms
For the algorithms with two parameters we applied the heuristics to choose the value of the second parameter as described in Table 2. These heuristics were designed to achieve the maximum simplification rate for a given Hausdorff error threshold . It was established that for the algorithms of Lang and "Perpendicular distance" the optimal value of is 0.25 and for > 0.25 the simplification rate does not increase (however, the execution time of these simplification algorithms increases leading to additional overhead).
Evaluation results. Using the polygons from the dataset MPEG 7 CE-Shape-1 we have measured the execution time of each suitable simplification algorithm in relation to the Hausdorff error threshold (see Fig. 8). These measurements can be considered as a computational overhead related to the optimization step of our skeletonizing algorithm. In order to compare the quality of the simplification algorithms, we have measured the dependence of the simplification rate on the error threshold value .  Fig. 9 shows that in general the largest polygon simplification for a given is achieved by the algorithms LA and ZS, which have approximately identical dependency curves. Slightly smaller simplification is accomplished using the algorithms PD, VW and PD, which show also almost undistinguishable behavior (except for the algorithm VW, which overperforms all other algorithms for small values of < 0.002). The lowest simplification rates are attained by OP and RW algorithms with nearly identical dependency curves. As we can observe from Fig. 9 and Fig. 10, the fastest algorithms OP and RW achieve the smallest simplification rate and, therefore, cannot guarantee the fastest performance of the skeletonization algorithm. Thus, we have measured the execution time of the proposed skeletonization algorithm using the simplified polygons. The results are shown in Fig. 8.
In order to take into account the overhead execution time of the simplification algorithms, we have measured the total execution time of the skeletonization algorithm including the respective simplification routines (see Fig. 11). Fig. 11 allows us to choose the fastest version of the optimization heuristics. However, we need to consider that the heuristic optimization step might affect the accuracy of the resulting skeleton. Therefore, we have also calculated the error of the optimized skeletonization algorithm. Such error is measured as the Hausdorff distance between the ground truth skeleton and the result of the optimized skeletonization algorithm (see Fig. 10). Fig. 11 shows that Ramer-Douglas-Peucker (DP) and Visvalingam-Whyatt (VW) algorithms allow us to speed-up our skeletonization method to the greatest extent. These two algorithms (only) overperform the optimization-free approach (NO) in case of small values of ≤ 0.001. The skeletonization based on Opheim and Reumann-Witkam algorithms exposes the smallest error among the other approaches (see Fig. 12). However, for < 0.002 these algorithms have a large computational overhead eliminating the effect of the optimization. Therefore, it is reasonable to apply them only for > 0.002. Note that the difference between skeletonizing errors decreasing as becomes smaller (see Fig. 12). We have computed 2-sample t-test to validate the hypothesis that algorithms DP and VW produce different average skeletonization errors. The hypothesis testing (see Table 4) showed that the skeletonizing errors produced by DP and VW are undistinguishable. Another hypothesis testing was performed to distinguish between the execution time between DP and VW algorithms. Table 5 and Fig. 11 show that for the most of the cases (except = 0.001) DP algorithm executes faster than VW.

DISCUSSION
Speed-accuracy trade-off. Since none of the tested algorithms minimizes the accuracy and execution time of the skeletonizing method at the same time, the optimal heuristic choice should base on the trade-off between accuracy and the speed. Thus, based on the performed computational experiments the following conclusions can be drawn: 1. If accuracy of the resulting skeleton is critical, then for > 0.002 the optimization can be performed using OP or RW algorithms. However, for < 0.002 the only reasonable optimization is using the DP or VW algorithms; 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]. Fig. 13 shows that for large values of simplification heuristics tends to regularize shape of the object in a way that the branches of the skeleton corresponding to small shape perturbation disappear. 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 brunches of skeleton.   = .
for the top row of images, = . for the bottom row of images

CONCLUSION
In this paper, we proposed an optimized algorithm for computing the Voronoi skeleton based on polygonal data. This topic is of relevance because of its direct relation to the optimization tasks in image processing and computer graphics (in particular, the efficient processing of vectorized images). We have illustrated in detail the main steps of the proposed skeletonization algorithm. It was established that the complexity of the algorithm is O (N log N), where N is the number of vertices in a polygon. We have also proposed theoretically justified optimization heuristic, which is based on polygon/polyline simplification algorithms. In order to evaluate and prove the efficiency of such heuristic, a series of computational experiments were conducted based on the polygons obtained from MPEG 7 CE-Shape-1 dataset, which represent the most commonly observed shapes in computer graphics and vision. In order to determine the most suitable optimization heuristic, we have evaluated seven different appropriate state-of-the-art simplification algorithms. We have measured the execution time of the skeletonization algorithm with and without the optimization and determined the computational overheads related to such optimizations. Also, we determined the accuracy of the optimized skeletonization algorithm depending on the applied optimization. As a result, we have established the criteria, which allow us to choose the optimal heuristics depending on the system's requirements.