{'fraction_non_alphanumeric': 0.0296875, 'fraction_numerical': 0.009114583333333334, 'mean_word_length': 4.332664712235766, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing a Compact Spline Representation of the Medial Axis Transform of a 2D Shape

We present a full pipeline for computing the medial axis transform of an arbitrary 2D shape. The instability of the medial axis transform is overcome by a pruning algorithm guided by a user-defined Hausdorff distance threshold. The stable medial axis transform is then approximated by spline curves in 3D to produce a smooth and compact representation. These spline curves are computed by minimizing the approximation error between the input shape and the shape represented by the medial axis transform. Our results on various 2D shapes suggest that our method is practical and effective, and yields faithful and compact representations of medial axis transforms of 2D shapes.


Introduction
The notion of the medial axis transform was first introduced by Blum [1] as an intrinsic shape representation. The medial axis of an object O is the set of interior points having at least two closest points on the boundary ∂O of O. In the 2D space, each point on the medial axis is the center of a circle, namely a medial circle, which is the maximal inscribed circle contained in O and tangent to ∂O in at least two points. To encode the complete shape information of the object, each point on the medial axis is assigned with the radius, which could be 0, of its associated medial circle. Therefore, a radius function could be defined on the medial axis. The medial axis coupled with a radius function is referred to as the medial axis transform (MAT). Each point in the MAT, called medial point, has three dimensions, which indicates its 2D position and the radius. The MAT is a complete shape representation in the sense that the object boundary can be reconstructed exactly from its MAT as the envelope of all the medial circles.
The MAT encodes rich information of a shape, such as local thickness, symmetry and its part structure, which is not possessed by alternative boundary surface representations. Therefore the MAT has been used extensively in a wide spectrum of applications, including shape analysis [2], shape deformation [3] and artistic rendering [4]. Detailed discussions on properties and applications of the MAT can be found in the book [5].
The MAT, on the other hand, is well-known suffering from the instability problem: small variations of the shape boundary may yield a large change to its MAT. While boundary noise is ubiquitous in data acquisition due to errors introduced in scanning, sampling and other numerical processing, the medial axis thus often has excessive geometric complexity and pathological topology, rendering it generally useless in practice unless it is cleaned up. A lot of existing algorithms have been developed to resolve the instability issue. As a common practice, unstable branches of the MAT induced by boundary noise are pruned based on certain measures [6,7,8]. Different criteria have been introduced to characterize the difference between the original shape and the reconstructed shape, or describe the remaining stable MAT with some intrinsic measures [9]. These methods focus on simplifying the topology of the MAT by pruning unstable branches, producing a topologically clean MAT which is nevertheless still represented by a large number of sample medial points. This is partially due to the prevailing choice taking the union of the sample medial circles (or medial spheres in 3D) as an approximation when evaluating the approximation error during both pruning stage and shape reconstruction stage. As a consequence, a large number of medial circles are often necessary to attain a good approximation to the MAT [10,11].
To achieve a smooth representation and further reduce the geometric complexity of the MAT, we propose to represent the MAT as spline curves. See Fig. 1 for a comparison of the two shape representations. The medial axis also possesses piecewise C 2 continuity in each medial branch. Although a similar smooth representation of the MAT has been used for modeling and segmentation purposes [12], a fully automatic way of obtaining a smooth MAT for an arbitrary shape is still missing.
In this work, we propose a complete pipeline which automatically computes a stable and compact medial axis  Figure 1: Given a shape swept by a moving circle, its medial axis transform is approximated as a sample set of 57 medial circles. (a) The approximate shape is the union of the medial circles. (b) On the other hand, to achieve the same accuracy with spline medial axis approximation, we need only 5 control points. The control polygon of the spline curve is shown in red.
transform which accurately approximates an arbitrary 2D shape. Given an error thresholdε, our algorithm guarantees that the Hausdorff distance between the boundary of the original shape and the boundary of the reconstructed shape is at mostε. Our method involves the pruning of unstable medial branches with an error-driven filtering process and the computation of a compact and accurate spline approximation to the MAT.
Compared to other works on medial axis computation, our algorithm possesses the following advantages: • Topological filtering with error control: Our filtering process is guided by a user-defined error thresholdε to ensure approximation accuracy while removing noisy, unstable branches as much as possible.
• Computing a compact geometric representation: We use spline curves to approximate the MAT, resulting in a highly compact representation. An optimization process is developed to make sure that the reconstructed boundary best fits the input shape, meeting a user specified error tolerance.
This paper is organized as follows. We start with a brief review of the previous work related to medial axis computation in Section 2. We then define the piecewise smooth medial representation in Section 3. In Section 4, we introduce our main algorithm; the implementation details are then provided in Section 5. We present experimental results in Section 6 and finally conclude the paper in Section 7.

Related Work
There is a vast amount of research studies about medial axis computation and representation. Here, we will review only those which are in close relation to our work.

Medial axis computation
Exact medial axis computation is possible only for simple or special shapes, such as polyhedra [13,14]. For free-form shapes, medial axis approximations are widely used in practice. There are several main approaches to computing the medial axis approximation: pixel or voxel-based methods that compute the medial axis using a thinning operation [15]; methods based on distance transform [16,17,18,19], often performed on a regular or adaptive grid; the divide-and-conquer methods [20], performed on spline curve boundaries; the tracing approaches [21], by tracing along the shape boundary or the seam curves; and the Voronoi diagram (VD) based methods [6,22,23,24,25].
Among these, the VD based approach stands out due to its theoretical guarantee and efficient computation. As a preprocessing step, we obtain an initial discrete medial axis of a shape using the VD based algorithm. The VD based method assumes that the boundary of an input shape O is a smooth curve and is sampled by a dense discrete set P of points (Fig. 2a), with the sampling density determined by the local feature size [22] in order to capture the boundary topology correctly. The Voronoi diagram of P is computed and the Voronoi vertices interior to O are taken to approximate the medial axis of Voronoi diagram (Fig. 2b), since a point on the Voronoi diagram is also characterized by having at least two closest points among the sample points.

Handling instability
Many studies have been conducted to understand and resolve the instability problem of the MAT. We review here several typical methods, whereas a survey can be found in [9]. One general approach is to define certain measures for the significance of a medial point, and to filter medial points against a user-defined threshold, thereby removing unstable branches of the medial axis. Examples include the angle-based methods which consider the separation angle or the object angle (i.e., the angle spanned by the closest contacting points) [16,26,27] and scaled axis transform (SAT) which essentially exploits the rate of change of the radius function as the filtering condition [7].
Another approach of computing a stable MAT is to consider the difference between the initial shape and the reconstructed shape from the pruned MAT [9]. The filtering The reconstructed shape from the medial axis transform can be obtained from our envelope reconstruction of a densely sampled medial axis and its associated medial circles.
step in our algorithm resembles this latter approach by considering the Hausdorff distance between the boundary of the input and the approximate shape to ensure the approximation accuracy of the output stable MAT.

Smooth medial axis representation
Yushkevich et al. [12] propose a continuous medial representation by modeling the MAT with cubic B-splines, as an extension of its discrete counterpart called the mrep [28]. The m-rep is built upon a sparse set of medial atoms, each of which encapsulates the 2D position of a medial point and the corresponding spoke vectors from the 2D medial point to the closest points on the object boundary. The continuous m-rep [12] (cm-rep), on the other hand, uses control points in cubic B-splines to describe the MAT, which must meet specific constraints defined on the implied boundary. In applying cm-reps to object modeling and image segmentation, a template cm-rep model is first built manually which is then deformed to fit a target shape. There is currently no method for automatically computing a smooth curve approximation to the MAT of a 2D shape.
The smooth medial representation we propose in this work is inspired by the cm-rep. However, instead of medial atoms, we adopt the simple medial points as our basic control entities. While the shape boundary is explicitly given by the trace of the spoke vectors for a cm-rep, the shape represented by our smooth medial axis transform is implicitly described by the union of the envelopes of adjacent medial circles. Using medial circles not only provides a more compact representation but also relieves the burden of ensuring the shape boundary as defined by the spoke vectors is intact and consistent, without compromising the accuracy in shape approximation.

Preliminaries
In this section we will introduce the notation for describing our method. Refer to Fig. 3 for an illustration of the medial axis transform of a 2D shape. Consider a 2D shape O as a compact connected subset of R 2 with boundary ∂O (Fig. 3a). The medial axis MA of O is defined as the set of points inside O with at least two nearest neighbors on ∂O. Medial circles could be described from points in MA, associated with the local thickness (radius) values. The medial axis transform MAT of O is the set of such medial circles. Each point in MAT , which we call a medial point, is specified by a 3D vector v i = (u i , r i ), where u i = (x i , y i ) ∈ R 2 denotes the medial position and r i ∈ R is the distance from u i to the shape boundary ∂O. In other words, v i describes a 2D medial circle, centered at u i with radius r i , which is tangent to ∂O in at least two points. In the extreme case, for example, a sharp corner of ∂O, the radius is actually 0.
A medial point v i tangent to ∂O at two distinct points is called a regular medial point. A branch is a maximal curve segment comprising regular medial points. A joint is a medial point at which three or more branches meet, and hence a joint has three or more tangent points with ∂O. An end-point is a medial point corresponding to an oscillating circle of ∂O.
Given a 2D shape O, we consider an approximation to its MAT by a piecewise smooth medial axis transform, denoted M = {C j } that consists of a set of medial spline curves connected at the joints of MAT . Each spline curve C j parameterizes both the position and the radius of medial points. Hence, any point on C j gives a medial point on M. Similarly, each C j is a maximal curve segment comprising regular medial points and it corresponds to a branch in MAT . An end-point of C j is either an end-point of M, or a joint of M at which three of more branches meet (Fig. 3b).
The envelope C j of C j is specified as the union of all medial circles on C j , noted that C j is a continuous curve. Let M be the union of { C j }. Then M represents a shape approximation to O. Compared with discrete medial axis representation in which a 2D shape is approximated using a union of medial circles sampled on MAT , the piecewise smooth medial axis transform provides a compact yet accurate shape representation. While a set of densely sampled medial circles needs to be sampled on a branch of MAT to well approximate the shape the branch represents ( Fig. 3c), it takes far fewer control points to define spline curves to accurately approximate the medial axis transform.

Computational Framework
In this section we shall present a full pipeline for automatically generating a piecewise smooth medial axis transform for an arbitrary 2D geometric shape O with smooth boundary ∂O. To start, we assume that O is sampled by a set of dense sample points {p i ∈ R 2 } on ∂O, meeting the sampling conditions [22] to faithfully capture the topology and features of the boundary ∂O.
Our algorithm consists of two main steps: (1) noise filtering for generating a stable medial axis transform M s , and (2) geometric simplification for computing a splinebased medial axis transform M based on the branch structure of M s . Our goal is to obtain M, whose reconstructed shape M well approximates the original shape O, and both the above steps are governed by a quantitative measure ε(O, M) of the approximation error of M to O.
While the differences of two shapes are commonly measured by a Hausdorff distance (by considering the shapes as point sets), we use instead the one-sided Hausdorff distance from the original boundary ∂O to the reconstructed boundary ∂ M of M as the error measure ε(O, M), which is easier to compute. As ∂O is generally more complex than the simplified ∂ M, the one-sided Hausdorff distance is usually a good approximation to the Hausdorff distance of the boundaries of the two shapes.
Our algorithm is outlined in Algorithm 1. We shall present the details in subsequent subsections.

Noise filtering
In this step, we aim at generating a stable medial axis transform M s for an input 2D shape O. Furthermore, we ensure that the approximation error of the reconstructed shape M s to O, denoted as ε(O, M s ), is no larger than a user specified error thresholdε. In our framework, ε(O, M s ) is measured by the one-sided Hausdorff distance from ∂O to the boundary ∂ M s of M s . By noise pruning, we shall obtain a stable medial axis transform M s that will be used as an input in the next step for computing M, a smooth spline representation of the medial axis transform.
Given a dense set of points {p i } sampled on the boundary ∂O of O, we first apply the Voronoi diagram based method [24] to obtain an initial medial axis transform The one-sided Hausdorff distance from ∂O to ∂ M s is defined as: where d(p i , ∂ M s ) is the Euclidean distance from boundary point p i to ∂ M s . The reconstructed shape M s is represented as the union of all medial circles in M s . Since the medial axis of M s is a subset of the Voronoi diagram of {p i }, p i always lies on or outside medial circle (u j , r j ). Hence, d(p i , ∂ M s ) equals the distance from p i to its nearest medial circle (u j , r j ):

Filtering strategy
To ensure that M s has the same homotopy as the initial medial axis M 0 , only the end-points of the medial axis transform might be pruned. The initial M 0 gives the exact reconstruction of the shape with respect to the sample points {p i } and hence ε(O, M 0 ) = 0. At each iteration of the filtering process, we check the error induced by the removal of a medial end-point v j , which equals If the induced error δ j ≥ε, v j will be treated as a feature of the medial axis and kept in M s . Otherwise, we treat v j as noise and prune it from M s . The process is repeated until all endpoints have been checked. Algorithm 2 describes the steps of our noise pruning algorithm.
The result of the above pruning procedure is in general not unique, since it might depend on the order of pruning. Pop an end-point v j from Q.

5:
If δ j <ε, remove v j , as well as its neighboring segment from M s and push the neighbor point of v j into Q if it becomes an end-point. 6: end while 7: return the stable medial axis transform M s . Nevertheless any such result will meet the specified error toleranceε, and can be used at the next stage of geometric simplification as an input to further compute a compact piecewise smooth medial axis transform.

Geometric simplification
After noise pruning, a stable medial axis transform is obtained. However, the pruning step removes only noisy medial points and leaves numerous discrete medial points on the stable medial branches. While those medial points are critical to a faithful representation for union of circles, in the envelope representation which we adopt, most of them contribute quite little to the shape and are thus nearly redundant. To get a concise medial axis transform, we should further reduce the number of medial points by utilizing a piecewise continuous representation. Such medial point decimation is driven by a user-specified approximation error thresholdε to ensure the approximation error is under control.
First, we construct an undirected graph G to represent the medial points filtered by noise pruning algorithm. Each medial point corresponds to one vertex in G. Two vertices in G are connected by an edge if their corresponding medial points are neighbors in the medial axis transform. We then obtain chains from the graph by grouping all edges which are connected without passing through a joint. Since each chain H j is single connected, we fit a smooth curve C j for it. By considering the smoothness, the ease of implementation and the representation ability, we choose the cubic B-spline to represent the smooth curve for each chain.

Piecewise cubic B-spline medial axis transformation initialization
We first sample dense points from the medial axis transform and find an initial cubic B-spline medial axis transform with open-uniform knot vectors to fit these medial points [29]. The reason that we choose this fitting algorithm is, it could generate good fitting curves only based on the number of control points; no initial curve is required. To ensure that the B-spline curve passes through the end-points, knots of multiplicity four are used at each end-point.
The fitting error of the C j curve is measured by the maximal Euclidean distance from medial points on chain H j to their projections on the fitting curve C j . Let f e be the maximal fitting error over all spline curves. Jiří Kosinka and Bert Jüttler [30] proved that, the boundary approximation error ε(O, M), measured by the one-sided Hausdorff distance from the original boundary ∂O to the reconstructed boundary ∂ M from M, is upper bounded by √ 2 · f e . Hence, we mark a spline curve C j as reliable, if the fitting error of C j does not exceedε/ √ 2. For every branch H j of M s , our fitting strategy is to search a reliable cubic B-spline curve with the minimal number of control points.

Medial axis transform optimization
Then, the distance between ∂O and ∂ M is minimized by optimizing the coordinates of B-spline control points. We'll define an energy function E(X) to formulate this problem. In the following presentation, the variables, where u i and r i denote the position and the radius of a control point.
The distance between the two curves is approximated by the sum of the squared distances.
where d 2 (p i , ∂ M) is the squared Euclidean distance from p i to ∂ M. Since M is the union of all C j , d 2 (p i , ∂ M) could be computed by the squared distance from p i to the boundary of its closest envelope C j , where ∂ C j represents the boundary of an envelope C j . As C j is the union of all medial circles it represents, a bottom-up way to compute d 2 (p i , ∂ C j ) is designed in our algorithm. By evenly and densely sampling in the spline parameter space, a set of piecewise linear segments are obtained. Let S be a segment sampled on C j . Then its reconstruction S is the envelope of its two medial circles, with a boundary ∂ S (see Fig. A.11). If the sampling density in spline parameter space is high enough, C j can be approximated accurately as a union of all such envelopes S. We consider the footpoint of p i on one envelope S as invalid, if the footpoint lies in the interior of another envelope. In that case, d 2 (p i , ∂ C j ) is the minimal valid d 2 (p i , ∂ S). We will provide the details for computing d 2 (p i , ∂ S) in Appendix A.
By minimizing E, an optimized M is obtained. Intuitively, more control points in M will lead to a smaller E and a smaller ε(O, M). As we know, the energy is defined in the L 2 space, while the approximation error ε(O, M) is defined in the L ∞ space. It is possible that even if E converges to a minimal value, ε(O, M) is still larger than ε. If that happens, our strategy is to insert new control points in M until ε(O, M) is no greater thanε.
In our algorithm, ε(O, M) is compared withε, at the end of the optimization process. If ε(O, M) >ε, our algorithm will pick up the boundary point p k with the largest distance to ∂ M, and insert a new control point to the spline curve whose envelope p k 's footpoint lies on. Another round of optimization will be executed, and these steps are repeated until a medial axis transform with a satisfactory approximation error is achieved. The steps are described in Algorithm 3. When inserting a new control point to a spline curve, the user is not required to specify the location of the new control point, but the spline curve fitting algorithm will compute and adjust the location based on the number of control points, which is similar to the spline curve initialization.

Implementation Details
In this section, some implementation details will be discussed.
Initial sampling. To get a good approximation of the medial axis transform, the input sample points {p i } should capture the boundary topology faithfully, as well as boundary features. The sampling condition we adopt is the local feature size condition discussed in [22]. The local feature size at a boundary point p i is the Euclidean distance from p i to its nearest point on the medial axis. In our case, if the boundary is r-sampled for r ≤ 0.25, we consider the set {p i } as a valid input.
Optimization. We apply the L-BFGS method [31], an iterative quasi-Newton method, to minimize the energy function E(X). The L-BFGS method takes the control points X in M as the input, evaluates the objective function E(X) and its gradient ∇E(X) in the optimization. Since no Hessian is involved in the computation, the L-BFGS method is quite efficient in finding the optimized medial axis transform.
Global check and local check. There're overlaps among envelopes. A footpoint is only valid if it is not contained in any other envelope. To determine the validity of a footpoint, a global check is required to traverse all envelopes, which is quite time-consuming as a repeated step in each iteration. We find that generally the footpoint will move only marginally after each iteration. The observation inspires us to adopt a hybrid check strategy, which uses the local check and calibrates the footpoint projection with a global check after a certain number of iterations.
Medial axis transform validation. The case where one medial circle is completely contained by another medial circle is not allowed in a valid medial axis transform.
Let v 1 = (x 1 , y 1 , r 1 ), v 2 = (x 2 , y 2 , r 2 ) be any two different medial points in a medial axis transform, and ∆x = x 1 − x 2 , ∆y = y 1 − y 2 , ∆r = |r 1 − r 2 |. In a valid medial axis transform, the slope of the two points, defined as tan α = ∆r should satisfy tan α ≤ 1. Although we do not incorporate this validity constraint in our optimization for the sake of simplicity, we have checked our computed examples and found no violation. However, to guarantee the validity explicitly, it is an interesting future problem to implement efficient constrained optimization to include this validity condition.

Experiments and Discussions
In this section, we choose some data sets and show how our algorithm performs on these shapes. All experiments are conducted on a Xeon 3.33 GHz PC with 12 GB RAM. We implement our system in C++ and adopt CGAL [32] to compute the Delaunay triangulation and generate the Voronoi diagram of sampling points on the shape boundary. The L ∞ approximation error ε and error threshold ε are normalized by scaling the diagonal of the bounding box of the input shape to 1.
In all figures presented in this paper, the blue regions are the reconstructed shapes from medial axis transforms and the black contours are boundaries of the input shapes. Black dots in spline medial axis transforms represent the positions of the B-spline control points. In other discrete medial axis transforms, black dots are the centers of medial circles.

Workflow
The workflow of our algorithm, from the initial medial axis transform to its piecewise smooth medial axis transform, is shown in Fig. 4. The initial medial axis transform is computed based on the Voronoi diagram of boundary points (Fig. 4a). Filtering the noise in Fig. 4a with an error thresholdε = 0.40%, we obtain a clean but discrete sample points of the medial axis transform (Fig. 4b), with 357 medial points and approximation error ε = 0.40%. In Fig. 4c,  12 cubic B-spline curves with 46 control points are used to fit medial points in the discrete medial axis transform, and an initial piecewise smooth medial axis transform is ready. After optimization, a piecewise cubic B-spline medial axis transform, whose approximation error ε = 0.27%, is shown in Fig. 4d. This example shows that the spline representation computed by our algorithm not only provides a smooth and compact representations for medial axis transforms, but can also improve the approximation accuracy.

Experiments on noise pruning
Our experiments demonstrate that the noise filtering algorithm (Algorithm 2) provides a reliable and stable medial axis transform with guaranteed approximation error. An examples is shown in Fig. 5. In Fig. 5b, the noisy branches in the initial medial axis have been filtered successfully. Meanwhile, the approximation error of the stable medial axis transform, 0.09%, is less than the specified error threshold 0.10%.
When handling shapes with perturbations on the boundary, our pruning algorithm also performs well. Although the initial medial axis in Fig. 6a contains numerous undesired branches, the resulting medial axis in Fig. 6b is topologically clean and stable. More noise pruning results are displayed in Fig. 10. This fact shows that our noise pruning algorithm is a convincing method to deal with noisy branches in the medial axis transform.
The error threshold can be fine-tuned to obtain different levels of details to be preserved. If the user would like to preserve more details of the input shape, a smaller error threshold should be applied, as shown in Fig. 7b. When a relatively large error threshold is applied, the main skeleton of the shape becomes clear, as shown in Fig. 7c.

Comparisons with other pruning strategies
There are existing medial axis pruning methods to filter noise from the initial medial axis transform. We perform a comparison of our method with two typical methods -the angle-based method and the scale axis transform (SAT) on the butterfly shape in Fig. 8. The medial axis of a butterfly with a smooth boundary (Fig. 8a) is chosen as the ground truth. Fig. 8b is generated after adding white noise on the smooth boundary in Fig. 8a, which results in a medial axis with many unstable branches. Fig. 8c to 8e show filtering results on Fig. 8b by angle-based method [8] with angle threshold 0.81, SAT [7] with scale parameter 1.12 and our method withε = 0.10%, respectively.
It can be observed that the result of our method is the most close to the ground truth. The angle-based method fails to prune the noisy branch in upper right wing, shown in the red rectangle in Fig. 8c. The SAT method ignores a small noisy branch in the tail, shown in the zoom-in view, while a branch of the right wing has been filtered erroneously as shown in Fig. 8d, highlighted with a red dash rectangle. Our algorithm keeps the skeleton properly, as shown in the green rectangle in Fig. 8e. Besides its power in pruning noisy branches, our result is also most accurate by achieving an approximation error 0.1%, much less than 0.34%, the error of the other two methods.

Advantages of spline representation
As mentioned before, spline curves provide a compact and smooth representation of medial axis transforms. Meanwhile, the control points of spline curves make the medial axis transform much easier to manipulate in applications such as shape editing and deformation. Fig. 9 shows comparisons between discrete medial axis transforms and the corresponding piecewise cubic B-spline medial axis transforms. The number of control points are largely reduced in the spline representation, compared with the number of medial points in the discrete medial axis transform, while the approximation error of the spline representation is comparable with that of the discrete medial axis transform.  Table 1: Comparisons of discrete medial axis transforms and spline medial axis transforms (shapes of Fig. 9). #pi is the number of boundary points. ε is the approximation error of medial reconstruction to O. Ms represents the discrete medial axis transform, and M is the corresponding spline medial axis transform. |V | is the number of medial/control points in Ms/M. Compactness records the medial points reduction from Ms to M.

Performances
We tested our algorithm on various 2D shapes. Some of the results are shown in Fig. 10. In each set of shapes, three medial axis transforms are displayed: the initial M 0 generated from a Voronoi-based approach; a stable medial axis transform M s computed by Algorithm 2 and a compact spline-based medial axis transform M optimized by Algorithm 3. These examples demonstrate that our algorithm provides an efficient way to obtain faithful and compact medial axis transforms for different 2D shapes, following the control of specified error thresholds. The statistical analysis in Table 2 shows that the number of points required to approximate the medial axis transform is greatly scaled down with our spline representation of the medial axis transform. In addition, the L-BFGS method assures the efficiency of our optimization process.

Conclusion
We propose a framework for computing an accurate and compact medial axis transform for an arbitrary 2D shape.
The approximation quality of the medial axis transform is guaranteed to be less than a user specified threshold. Moreover, users have the freedom to choose different types of splines to represent the medial axis transform according to their requirements and applications. A noise pruning algorithm for the medial axis transform is also proposed and integrated in the framework. This noise pruning algorithm filters noise in the medial axis transform robustly, and provides a good initial medial axis transform for optimization in our framework. Experimental results confirm the effectiveness of our algorithm. In future, we plan to study how to apply the new spline-based representation of the medial axis transform in shape deformation and shape matching, and extend the framework for the computation of the medial axis transform of a 3D volume.

Appendix A. Error term computation
As discussed in Section 4.2, the energy function E(X) (Equation 3) could be computed from d 2 (p i , ∂ S), which is the squared Euclidean distance from boundary point p i to the envelope boundary of a segment S sampled on the curves of M. For each segment S, there are two points v 1 and v 2 , which represent two medial circles (u 1 , r 1 ) and (u 1 , r 1 ) in 2D space. And there are one or two common external tangent for two medial circles, which is implied by the property of the medial axis transform.
Consider the case where two circles have two distinct external tangent lines (Fig. A.11). Let q i (i = 1, 2, 3, 4) be the tangent points. The external tangent lines through q i are represented as L 12 and L 34 . The circle is partitioned into two circular arcs by tangent points, and the two outer arcs are denoted as A 14 and A 23 , respectively. The boundary of the shape represented by S, denoted by ∂ S, consists of two outer arcs A 14 , A 23 and two external tangent lines     Table 2: Performances of our algorithm (shapes of Fig. 10). #pi is the number of boundary points.ε is the user-definded error threshold. ε is the approximation error of medial reconstruction to O. |V | is the number of medial/control points in medial axes. M0 is the initial medial axis transform based on a Voronoi-based approach. Ms represents the stable medial axis transform computed by our noise pruning algorithm, and M is the spline medial axis transform optimized by our geometric simplification algorithm. Compactness shows the medial points reduction from Ms to M. Time column records the time cost of geometric simplification from Ms to M. L 12 , L 34 . Let β 1 = ∠p i u 1 u 2 , β 2 = ∠p i u 2 u 1 , α 1 represents ∠q 1 u 1 u 2 or ∠q 4 u 1 u 2 , α 2 be ∠q 2 u 2 u 1 or ∠q 3 u 2 u 1 . By comparing these angles, we can easily tell which part p i 's footpoint lies in. In sequence, the distance from p i to its footpoint on ∂ S is computed as if β 1 < α 1 and β 2 < α 2 and p i , u 1 , u 2 are in clockwise order d 2 (p i , L 12 ) otherwise , where d 2 (p i , L 12 ) and d 2 (p i , L 34 ) are easy to compute, since we know the algebraic representation of two outer tangents, and the gradient of d 2 (p i , ∂ S) can be obtained similarly.
It is also possible that only one common external tangent exists for two medial circles, which occurs when the slope of two medial points equals 1, computed as Equation 5. In that case, S is a circle and computation of d 2 (p i , ∂ S) is trivial.