- Email: [email protected]
Medial axis of a planar region by offset self-intersections
Computer-Aided Design 41 (2009) 1050–1059
Contents lists available at ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Medial axis of a planar region by offset self-intersections R. Dorado ∗ Department of Mechanical and Mining Engineering, University of Jaén, Edificio A-3, Campus Las Lagunillas, 23071 Jaén, Spain
article
info
Article history: Received 31 December 2008 Accepted 22 August 2009 Keywords: Medial axis Medial axis transform Offset self-intersection Interpolation
abstract The medial axis (MA) of a planar region is the locus of those maximal disks contained within its boundary. This entity has many CAD/CAM applications. Approximations based on the Voronoi diagram are efficient for linear-arc boundaries, but such constructions are more difficult if the boundary is free. This paper proposes an algorithm for free-form boundaries that uses the relation between MA and offsets. It takes the curvature information from the boundary in order to find the self-intersections of successive offset curves. These self-intersection points belong to the MA and can be interpolated to obtain an approximation in Bézier form. This method also approximates the medial axis transform by using the offset distance to each self-intersection. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction
a The Medial Axis (MA) or skeleton of a planar region is the locus of centers of maximal disks contained within its boundary. Together with the function that gives the radius of each disk, the MA defines a Medial Axis Transform (MAT). This concept, from the differential geometry, was proposed to describe biological shapes by Blum [1, 2], who also described two definition models, a maximal disk model looking for a circle that touches two different boundary points, and a grass–fire model taking the region’s boundary as initial fire and propagates within the region until it intersects itself. The MAT properties are very well-known because of the numerous studies about this concept, mainly for planar regions. Ramanathan and Gurumoorthy [3] enumerated these properties. The MAT simplifies regional representation, and retains the original information as Patrikalakis and Maekawa [4] recalled. This is interesting in many computer animation applications because changing MA or the radius function is enough to deform the shape of a region, to construct finite element meshes from MA [5], or during numerical control machining when defining an optimum tool path for operations such as pocketing [6]. Despite the numerous efficient methods to construct the MA of domains with linear-arc boundaries, the MA of a free-form region is still a problem [7]. One option is to approximate the boundary using straight lines, and compute the MA by one of the existing methods for polygons. The drawback is the topological error of the solution because MA is sensitive to boundary continuity (Fig. 1). This error is a problem in applications where the planar region has to be reconstructed from its MA, e.g. to deform a region via the MA. Studies of related topics, such as the Voronoi diagram and the bisector, are other ways of computing this problem. The
∗
Corresponding author. Tel.: +34 953 212439; fax: +34 953 212870. E-mail address: [email protected].
0010-4485/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2009.08.005
b
Fig. 1. MA computed from: (a) exact boundary (b) polygonal approximation.
MA of a free-boundary region can be computed using the Voronoi diagram [4,8], a dual of Delaunay triangulation [9]. Some scholars have studied the bisector to compute Voronoi diagrams. Held [10] for curvilinear polygons and Chou [11], Farouki and Ramamurthy [12], Ramamurthy and Farouki [13,14], and Seong, Cohen and Elber [15] for curved boundaries have obtained the Voronoi diagram by trimming bisectors. The previous methods use distance functions, as Seong et al. [15] that compute the bisector curve as an intersection of two distance surfaces, it is similar to the relation between untrimmed MA segments and the intersection between developable surfaces proposed by Pottman and Waller [16]. On the other hand, the trimming operation increases the process complexity. Moreover, if the boundary is free, neither the Voronoi diagram nor the MA is, in general, a subset of the other [13]. To obtain the MA you have to trim and/or add segments of the Voronoi diagram. Recently, scholars have been searching for other solutions. Choi, Choi and Moon [17] proposed a decomposition lemma to divide a general region into fundamental domains, where it is easy to compute the MA. They presented a computation method [18] based on the following idea; find the MA vertices, divide the domain, and compute normal MA points between two vertices with a method analogous to the maximal disk model.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
1051
b
a
Fig. 2. (a) MAT definition. (b) Criterion of minimum distance.
Ramathan and Gurumoorthy [3] used a tracing algorithm based on the evaluation of both distance and curvature criteria. Another tracing algorithm was proposed by Degen [19], who looked for the centers of maximal disks inside a planar region, and used the curvature information and the theory of dupin cyclides for his predictor/corrector algorithm. Based on the moving Frenet frames and Cesareo’s approach, Cao and Liu [7] discovered an interesting relation between two adjoining curves and the MA, and created a suitable tracing algorithm for solving the problem for free boundaries. Existing methods to compute the MA of a free-form region have, in general, a relatively low accuracy [7]. Moreover, most of them use a polygon approximation of the boundary, so the correction of topological errors is difficult to guarantee. An additional processing is also required to eliminate the artificial segments that may arise due to the discretization [3]. In this paper, the author proposes a method to approximate the MA of a free-boundary planar region. This method uses the exact analytic boundary, thus is free of the topologically incorrect segments produced by polygonal or curvilinear approximations. The solutions are MA segments and therefore, additional trimming operations, as in Voronoi techniques, are not needed. On the other hand, the approximation accuracy depends on the polynomial degree used to interpolate the computed MA points. To check the accuracy, a measure of the interpolation error is defined. The algorithm is analogous to the grass–fire model, because it finds the self-intersections of the successive boundary offset belonging to the MA. Note that grass–fire methods, like the proposed by Siddiqi et al. [20] or the methods reviewed by Damon [21], propagate a front through all pixels within a region of a digital image in order to choose possible candidates of MA, Oliveira and De Figueredo [22] developed a similar idea via a function distance. Consequently, it is difficult to guarantee that the solution is a thin set, i.e. without an interior, or is a connected set, and to compute the MA points accurately. By contrast, the proposed method uses a continuous definition of the boundary and its successive offsets instead of a digital image of the region, hence it guarantees an accurate, thin and connected (by interpolation) solution. Self-Intersection points are offset points with two or more different parameter values. For a parameter value u+ , it is possible to find a different parameter u− and the offset distance d for these parameter values, by solving a non-linear equation system. The solution for this system can be obtained using the Newton–Rapshon method with initial values of d and u− . These initial values can be approximated using information on the boundary curvature. In order to improve the algorithm and to deal with the least possible computations, this method interpolates the solutions at Chebyshev points. The solution is simple, natural, and thanks to the polynomial interpolation, with the same boundary
Fig. 3. Curvature sign and radius of curvature.
parameter. Finally, the algorithm gives the MAT, by interpolation of the offset self-intersection distance. This approximation is defined for regions with a general piecewise boundary of at least G1 (equal tangent direction at the piecewise unions). The extension to C 0 boundaries (equal points at the piecewise unions) is not the goal of this method, but Section 5.3 describes a simple example of this extension. The paper is arranged as follows. Section 2 describes the MA and offset relations. Section 3 shows a MA point classification and a definition of fundamental domains. It is possible to find the MA of general boundaries by the decomposition lemma [17]. If you can solve the fundamental case, you can solve the general one. Sections 4 and 5 describe the procedure for computing the MA points by offset self-intersections. The extension to general boundaries is based on the previous method and the minimum distance criterion, as shown in Section 6. Finally, some general examples and main conclusions are drawn in Sections 7 and 8 respectively. 2. Medial axis and offset self-intersections The MA or skeleton cs (u) of a planar region Ω can be expressed as a function of its boundary c(u), a radius function r (u) that gives the radius of maximal disks, and the unit normal directions from the boundary n(u). Let c(u), u ∈ [u0 , un ] be a parameterized closed boundary spline curve without intersections, in such a way that the region Ω is always to the left of parameter’s increasing direction, then the MA is given by (Fig. 2a): cs (u) = c(u) + r (u)n(u).
(1)
The unit normal n(u) and the curvature k(u) of c(u) are given by [4]:
(cu × cuu ) · ez , (2) |cu |3 where t(u) = cu / |cu | is the unit tangent of c(u), and ez = {0, 0, 1} n(u) = ez × t(u),
k(u) =
is a unit normal for the region plane. The positive curvature corresponds to c(u) points where the normal vector points out toward inside of the domain (Fig. 3).
1052
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
On the other hand, the MAT is the ordered set of all pair MA points and radius functions {cs (u), r (u)}. Note that MA and MAT depend on the radius function r (u). The Eq. (1) defines the centers of maximal disks. A disk is maximal when it is inside a region and does not contain other disks. According to Fig. 2b, the minimum distance dmin from the boundary to a disk center inside it and the disk radius R determine if the region contains the disk. These situations are possible and defining if a disk center belongs to the MA:
• dmin = R. Disks are within the region and have contact components (points or arcs) with the boundary. There are two cases: (i) One contact component. If the disk is an osculating one with the same tangent and curvature that the boundary (Fig. 3), its center is a MA point, because it is the maximal disk within the region with that contact (disk D1 in Fig. 2b). Otherwise, the disk is not maximal, so its center is not a MA point (disk D2 in Fig. 2b). (ii) More than one contact component. The centers are MA points, because they are the maximal disks within the region with those contact components (disk D3 in Fig. 2b). • dmin > R. This is the condition for all disks within the region without contacts with the boundary. They are not maximal, so their centers are not MA points (disk D4 in Fig. 2.b). • dmin < R. Disks intersect the boundary. They are not contained in the region and thus their centers are not MA points (disk D5 in Fig. 2b). The radius function r (u) defines the radii of maximal disks, and it is related to the offset distances. The proposed method uses this relation to compute the MA. In view of the applications for offset curves as tool path generation, feature recognition, or definition of tolerance regions, the previous idea is known, but it is used in the inverse problem to compute trimmed offsets. Chiang, Hoffmann and Lynch [23] use the MA to trim offsets, their idea is based on the fact that the MA cuts the offset of a closed curve at the self-intersections or at irregular offset points (isolated points and cusps). Choi et al. [24] show a completed study about the MA and offset relation. They note the coincidence between offset self-intersections and MA points, and use this relation to obtain a domain decomposition for offsets. Finally, Cao and Liu [7] also use the same idea to compute trimmed offset curves. The following definitions are useful to understand the relation between offset and MA. The offset cd (u) of a closed curve c(u) (with respect to distance d from the boundary to the interior of the region) is defined by: cd (u) = c(u) + dn(u).
(3)
The unit offset tangent td (u) and the offset curvature kd (u) depend on the tangent t(u) and the radius of curvature (radius of the osculating circle, Fig. 3) ρ(u) = 1/ |k(u)| of c(u): td ( u) = kd (u) =
1 − (d/ρ(u)) t(u), |1 − (d/ρ(u))|
(4)
1/ρ(u) . |1 − (d/ρ(u))|
(5)
Fig. 4. Relation between offset self-intersections and MA points.
Fig. 5. Decomposition lemma.
words, a self-intersection point is at the same distance respect to various boundary points and thus it belongs to the MA. Recall that the previous idea is used to trim offset curves in practical applications. This paper proposes the inverse problem, by computing the MA from offset self-intersections at different distances, as Fig. 4 shows. Instead of propagating the offset into the region, the algorithm solves a non-linear system of 2 equations for increasing parameter values u+ i , i = 0, 1, . . . , n. The solutions are the offset distances di to the MA and the parameter values u− i where the MA disks have a second contact point with the boundary. 3. Fundamental domains of a region Choi et al. [17,18] show that for a domain Ω with finite boundary curves linked with at least C 0 continuity (a spline curve), the MA is a connected geometric graph with a finite number of vertices and edges. Their decomposition lemma proves that a domain Ω can be divided into subdomains Ωj ,j = 0, 1, . . . , m, whose MAs are the Ω MA edges and vertices. A classification of MA points helps to obtain domain decomposition. According to Fig. 5, the MA points can be:
• Normal points. Center of maximal disks, of which each disk has two contact points with the boundary.
• Terminal points. These points are origin or end of one MA edge.
When distance d is equal to the radius of curvature, the offset has singularities (4) and (5). If d is larger than the minimum radius of curvature, the offset curve shows self-intersections. The Fig. 4 shows an example, note that ordinary MA points are centers of maximal disks with two contact points. The symbols +, − are used to distinguish between the parametric values of each point. Given the medial axis Eq. (1) and the offset Eq. (3), a self-intersection point with parameter values u+ , u− and offset distance d = r (u+ ) = r (u− ), lies on the medial axis. In other
(i) End-points are the centers of disks that have 1 contact component with the boundary: 1 point or 1 common circular arc. (ii) Corner points, boundary points where the unit tangent vector field is discontinuous. • Branch points. These points are the origin of more than 1 MA edges. (i) Non-generic points have 2 contact components with the boundary. (ii) Bifurcation points have 3 or more contact components with the boundary.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
1053
+ + boundary parameters u+ i , i = 0, 1, . . . , n, where u0 and un are the parameter values of the end-points, the corresponding di and u− i values to each MA point are the solution to the next selfintersection equations (+ and − are used to distinguish between the two contact components of the MA points): − cdi (u+ i ) = cdi (ui ),
(6)
the unknowns are di , ui . The non-linear system (6) can be solved by the Newton method, which has a quadratic convergence near a root (computational accuracy used < 10−8 ). The solutions u− i and + di for each u+ i give the MA points cs (ui ) by MA equation (1). Note − that to recover the boundary from the MA, u+ i and ui are necessary. The end-point parameters and radii of curvature give the initial values for solving the Eq. (6). The offset tangent and curvature equations (4) and (5) show that self-intersections begin when the offset distance is equal to the minimum radius of curvature ρ(u), thus self-intersection points start at the end-points computed in − Section 4.1. The MA computation begins with u+ 0 , u0 (for end+ + − points with 1 contact point u = u ), d0 = ρ(u0 ) = r (u+ 0 ) corresponding to one end-point. For u+ and the initial values i + − u− i−1 , di−1 = ρ(ui−1 ), the solutions for Eq. (6) are ui , di . Repeat this process until the parameter value corresponds to the other end+ + point u+ n , dn = ρ(un ) = r (un ). Finally, polynomial interpolation of MA points cs (u+ i ) and offset distances di = r (u+ ) gives the approximations of MA i and the radius function respectively in Bézier form, the standard use in practical applications [25]. This representation is easy to incorporate into CAD software, is easier to implement than the rational one, and admits a rational representation with unit weights. The main problem for obtaining an efficient interpolation is how to choose the u+ i values. A polynomial interpolation of equally spaced points (u+ i ) is highly ill-conditioned: small changes in the data may cause huge changes in the interpolant [26]. A numerically stable approach is to use point-sets that are clustered at the end-points of the interval with an asymptotic density −
Fig. 6. End-points computation.
In the following sections, the term vertices is used instead of terminal or/and branch points when the context refers to both types of MA point. A fundamental domain is a subdomain whose skeleton only has 2 vertices (terminal or/and branch points), and normal points. They are the smallest subdomains of a region and, thus, suitable cases for showing the performance of a computation method for normal MA points. 4. MAT computation for fundamental domains This section describes how to get the terminal points of a fundamental domain and how to use them to get its MA edge. For a C 0 boundary the terminal points are the boundary convex corners. Otherwise, a fundamental domain has end-points. 4.1. End-points computation Fundamental domains have no branch points. If the boundary is a parametric spline curve of at least G1 , the MA is a curve with two end-points. These points correspond to the disks with radii of curvature minimum within the region and therefore as different research works note [11,18], the boundary curvature allows for computing them. From the normal and curvature definition (2), and the minimum distance criterion, the end-points are the center of the osculating circles where the positive curvature k(u) attains its maxima. To obtain the end-points: 1. Compute the parameters’ values u where the curvature attains its maxima. If the curvature has a discontinuity at u, choose the u value of the greater curvature. Fig. 6 shows the parameters selected for 2 different end-points. 2. Compute the minimum curvature radii ρ(u) as the inverse of maximal curvatures. 3. Evaluate whether the disks with these minimum radii are within the domain. The minimum distance from the boundary to the center of the disk must be equal to its radius. As Cao and Liu [7] propose, the minimum distance parameter belongs to a set of parameters where the tangent cu (u) is perpendicular to the direction between the disk center p and the boundary c(u). The minimum distance associated with the roots of (p − c(u)) · cu (u) = 0 is the minimum distance needed. 4. The selected centers of the disks within the domain are the endpoints. Note that non-generic points are also MA points with maximum curvature. The number of contact components indicates if a point computed by the previous algorithm is non-generic. 4.2. MAT computation by offset self-intersections A self-intersection point is also a point of the MA, since the offset distance d, the parameter value u at this point, and the radius function are related r (u) = d. For a set of ordered
−1/2
proportional to 1 − u2 as n + 1 → ∞. The families of Chebyshev points are examples of these clustered point sets. These points are obtained by projecting equally-spaced points on the unit circle down to the unit interval [−1, 1]. Therefore, to compute a well-conditioned interpolation in a general domain u ∈ [u0 , un ], the u+ i can be the set of points: +
ui = 1 − cos
iπ n
+ u+ n − u0
2
+ u+ 0 ,
i = 0, 1, . . . , n.
(7)
The previous equation is a linear transformation of the Chebyshev points. Therefore, the point set (7) has the same density, between the parameters of two end-points, that the Chebyshev points over the unit interval. The exact MAT solution (MA and the radius function) is, in general, unknown, so the interpolation accuracyis also unknown. On the other hand, a valid interpolation MAT set c∗s (u), r ∗ (u) has
to verify Eq. (1): n(u) · c∗s (u) − c(u) /r ∗ (u) = 1, so a way of measuring the interpolation error s is:
n(u) · c∗s (u) − c(u) s = 1 − . r ∗ (u)
(8)
If s is larger than a desired value, we need more interpolation data points. A similar procedure is used to recover the boundary by substituting the MAT interpolation at Eq. (1). The boundary approximation error c (u) can be measured by the distance between the
1054
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
b
a
Fig. 7. MAT of an ellipse region (a) Degree 2 approximation (b) Degree 5 approximation.
Fig. 8. MAT approximation of a fundamental domain bounded by a degree-3 B-spline curve. (a) Degree 4 approximation (b) Degree 6 approximation.
boundary c(u) and the approximation c∗ (u). These curves have a common parameterization because of the approximation method.
c (u) = c(u) − c∗ (u) .
(9)
Steps to compute the MAT by offset self-intersections. 1. Compute the boundary end-points by the curvature information. For a general end-point with one contact component d0 = − + − ρ(u+ 0 ) = ρ(u0 ), and dn = ρ(un ) = ρ(un ). + 2. Choose the parameter values ui , i = 0, . . . , n with u+ 0 the parameter of one end-point and u+ n the other end-point parameter. Select a Chebyshev sequence (7).
3. Compute self-intersection points for each u+ i solving the Eq. (6). To compute the first self-intersection parameter u− 1 and d1 , + − − + + take u+ 1 and the initial values d1 = ρ(u0 ), u1 = u0 −(u1 − u0 ). The initial values for each new step are the solution of the + previous one until un . 4. Interpolate u+ , d i data points to obtain a Bézier MA and MAT i approximations. Finally, recover the boundary with the MAT set. In the following examples, the accuracy used to compute the curvature information and the self-intersections is 10−8 . The Figs. 7–14 were executed on a Mac with an Intel Core 2 Duo 2.8 GHz processor and 4GB RAM.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059 Table 1 CPU time used to compute Figs. 7–11.
a
1055
b
CPU time (s) Figure
End-points
MA edge
Total
7a 7b 8a 8b 9 10 11a 11b
0.018
0.008 0.020 0.025 0.061 0.010 0.030 0.010 0.070
0.026 0.038 0.055 0.091 0.220 0.290 0.030 0.270
0.030 0.210 0.260 0.020 0.200
Fig. 11. Degree 5 MAT approximation between different terminal points.
Fig. 9. Degree 2 MAT approximation between end-points with 1 contact arc. The boundary is a G1 spline curve.
End-points with 1 contact point have one parameter value. Therefore, the previous algorithm takes as initial values u+ 0 = u0 and u− = u . 0 0 A simple example is the MA of an ellipse, as shown in Fig. 7. Note that this method recaptures the exact solution for MA. On the other hand, the MAT and the recovered boundary are approximations. The previous result improves with more interpolation points. A general case is drawn at Fig. 8. The exact solution is unknown, so the interpolation accuracy can be measured by the interpolation error (8). The computational time is less than 0.1 s and depends on the interpolation degree (Table 1). The CPU time is mostly spent to compute the end-points. On the other hand, the number of interpolation points (degree interpolation + 1) increases the time to approximate the MA. 5.2. End-points with one contact circular arc A circular arc is common to those disks with center at these endpoints. The end-points are the MA’s for each circular arc. The initial − parameter values u+ 0 and u0 for computing the MA edge, are the arc terminal points. Fig. 9 is a simple example that shows how the offset collapses in the MA. There is a parallel pair of boundary segments in Fig. 9, note that Newton method gives one solution for each parameter value of the Chebyshev sequence along the boundary. Although the results are, in general, approximations, in this example the method provides an exact solution for MAT. The algorithm spends 0.22 s (Table 1) to compute a solution free of topological errors. 5.3. Domain with corners If a boundary region has convex corner points (with curvature
≥ 0 around them), the MA begins and/or ends at these points.
Fig. 10. Degree 5 MAT approximation for a boundary with convex corner points.
5. Examples of fundamental domains
For fundamental domains with 2 convex corner points, the initial − values for computing the MA are u+ 0 = u0 and u0 = u0 , where u0 is the parameter value of one corner. Fig. 10 shows an example. Note that concave corners are never on the MA, once they have been detected and isolated, the proposed algorithm can also be used in this case. Other examples of MA between different types of terminal points are drawn in Fig. 11. The CPU time is less than 0.3 s, and as the previous cases, the end-points consume almost all this time.
5.1. End-points with one contact point 6. Extension to general domains The previous method gives an approximation for different planar domains without branch points. The next examples show practical results and the extension to boundaries with corner points. Table 1 indicates the computational time.
Once the vertices (terminal or/and branch points) are known, the MA approximation of a general domain is the union of all fundamental MA computed between 2 successive vertices. The main
1056
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
Fig. 12. Degree 5 MAT approximation of a general region with 1 bifurcation point.
problem with general domains is computing the branch points, because the algorithm for fundamental domains helps to compute the end-points. Consider that we have two vertices and between them there is 1 unknown branch-point. Using the self-intersection algorithm between the first 2 vertices, some of the solutions are not MA points, because the minimum distance from the boundary to these points is less than the associated offset distance. The last correct MA point approximates the desired branch point. The previous idea helps to build an iteration procedure based on: firstly, compute an approximation between the 2 vertices and use the last correct MA point and the first wrong MA point to compute the next step, repeat the process until the difference between the 2 points computed in the last step is less than the desired precision, then take the last correct MA point as the branch
point. The results of this procedure are branch point coordinates and 2 parameter values where the associated disk has contact points with the boundary. The rest of the branch-point parameter values are where the boundary has the minimum distance with the branch-point. These steps summarize a self-intersecting algorithm for general domains that follows an analogous strategy to the Cao and Liu [7] one. 1. Define an ordered vector v with the couples of end and nongeneric parameters computed with the boundary curvature, as explained in Section 4.1. Sort the couples of parameters into smaller to greater values. The branch-point parameters are divided and ordered in couples from the first contact point until the same point again.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
1057
Table 2 CPU time used to compute the examples of Section 7. Figure
s <
CPU time (s) Vertices
12 13 14a
14b
14c
Fig. 13. Degree 20 MAT approximation of a general region with 2 bifurcation points.
2. Compute the self-intersections between the first 2 elements of v. Apply the minimum distance criterion. If their points do not belong to MA go to step 3. If all points belong to MA then go to step 4. 3. Compute the branch-point between the 2 first elements of v. Add the branch points’ parameters to v and compute the MA points between the 2 new first points of v. 4. Delete the 2 first elements of v and go to step 2. 5. Repeat from steps 2 to 4 until vector v is empty. 7. Examples of general domains The next examples help to understand the performance of the previous general algorithm. Fig. 12 shows the case of a domain with 3 end-points and 1 bifurcation point. To compute the MAT of this example:
• Generate a vector v with the end-points parameter values. = u− There are 3 end-points at u+ i , i = 0, 1, 2, so v = i + + + u0 , u1 , u2 . • Compute self-intersections between the 2 first elements of v. Compute the bifurcation point between the first 2 endpoints, as explained in Section 6. The bifurcation point has 3 contact components with the boundary ub0 , ub1 , ub2 . Join and sort the parameters to obtain a
+ + b b b new vector v = u+ 0 , u0 , u1 , u1 , u2 , u2 . In this example all MA segments can be computed from end-points with one contact, but if the construction begins in a bifurcation point, you need 2 initial parameters. • Compute the MA points by the offset self-intersections between the elements v and, finally, remove elements 1, 2 to obtain 1, 2 of + b b . a new v = u+ , u , u , u 1 1 2 2 There is no bifurcation point between the elements 1, 2 of v. • Compute the MA points between the elements 1, 2 of v. Next, b delete these elements to obtain a new v = u+ 2 , u2 . There is no bifurcation point between the elements 1, 2 of v, so repeat this step. • Compute the MA points between 1, 2 elements of v. Next, eliminate these elements to obtain a new v = {}. The vector v is empty.
A 2-branch point example obtained following the same general algorithm is drawn at Fig. 13. Finally, more difficult examples are drawn at Fig. 14. The procedure is the same for regions with holes inside (Fig. 14b), the main different is that the minimum distance criterion and the algorithm to compute the end-points also have to consider the holes. Fig. 14c shows the MA of a leaf computed from a photograph. Table 2 shows the CPU time used to compute the Figs. 12– 14. Note that branch and terminal points can be used to approximate the MA edges with different interpolation degrees and consequently their computational time is fixed. However, the
−3
11 × 10 8 × 10−7 4 × 10−2 6 × 10−4 2 × 10−7 10−2 2 × 10−4 5 × 10−7 10−2 10−4 3 × 10−7
1.00 1.90 3.71
4.41
20.80
MA edges
Total
0.05 0.81 0.17 0.41 2.52 0.22 0.51 3.83 0.84 1.99 14.61
1.05 2.71 3.88 4.12 6.23 4.63 4.92 8.24 21.64 22.79 35.41
number of interpolation points increases the CPU time, and the accuracy. The results are obtained in seconds for complex boundaries as in Fig. 14. The recent algorithm presented by Cao and Liu [7] spends minutes to compute the MA of regions with a similar number of end-branch points and MA edges. 8. Conclusions The method proposed approximates the MA of a free-form region using the boundary information. The algorithm has a general purpose. It can be of interest in tool-path optimization, and deformation of planar regions, because it is easy to program, easy to incorporate into CAD/CAM software, and it spends seconds to compute the examples shown in Sections 5 and 7. The MA points are computed by offset self-intersections. The relation between offset and MA is known, but it is, in general, used to trim offset curves [23,24] instead of computing the MA as in this work. Once the MA vertices (terminal or/and branch points) are known by the boundary curvature and the minimum distance criterion, their parameter and radius function provide the initial information needed to begin a tracing algorithm through the domain boundary parameters. This tracing algorithm is based on solving a non-linear system of two equations where the unknowns are the radius function and the second parameter associated with the computed MA point. Nonlinear systems are not a great drawback, because we know initial values near the solution. In order to obtain a more efficient approximation, the computed points are interpolated. The parameter values for the computed MA are the Chebyshev points, because they make the interpolation more stable. The solutions are curves in Bézier form being the standard in CAD/CAM, and have the same parameter as the boundary. The self-intersection method also gives the MAT of a region by interpolation of the radius function, and allows for recovering the boundary curve. The recovered boundary error is the distance in respect of the original curve, and the MA equation (1) allows for measuring the error of the MAT interpolation. The computational time depends on the boundary curve, because it determines the number of MA edges, vertices, and the interpolation degree needed to get accurate results. The CPU time used to compute the MA is less than 1 min in the examples (Table 2). Some differences between the proposed method and main existing constructions:
• Grass–fire methods are applied to digital images and thus it is difficult to guarantee thin and connected sets, and to compute accurate solutions. The proposed method provides an accurate and polynomial solution.
1058
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
b
a
c
Fig. 14. MAT of general regions, the boundary is a degree-3 B-spline (a) Hand (b) Hand with a hole (c) Oak leaf.
• Voronoi methods need additional trimming operations because if the boundary is free the Voronoi diagram is, in general, not the MA. Offset self-intersections are directly the MA points. • Maximal disks methods uses the differential properties of the MA and the boundary. The accuracy depends on the initial values, the tracing step size and the iterative step size. A corrector or predictor algorithm [19,7] improve the accuracy. The proposed algorithm is based on the relation between the offset distance and the MAT radius function. The accuracy depends on the computation of vertices and MA points (computational accuracy < 10−8 ) and the interpolation degree. Additional algorithms to improve the accuracy are not needed. Acknowledgements This research is supported by the Spanish Ministerio de Innovación y Ciencia, under research grants JC2008-00107, DPI200609091, DPI2009-10078, co-financed by the ERDF (European Regional Development Fund). The author thanks Dr. Borut Žalik for several helpful discussions on medial axis. References [1] Blum H. A transformation for extracting new descriptors of shape. In: Models for perception of speech and visual form. Weinant Wathen-Dunn MIT Press; 1967.
[2] Blum H. Biological shape and visual science (Part I). Journal of Theoretical Biology 1973;38:205–87. [3] Ramanathan M, Gurumoorthy B. Constructing medial axis transform of planar domains with curved boundaries. Computer-Aided Design 2002;35: 619–632. [4] Patrikalakis NM, Maekawa T. Shape interrogation for computer aided design and manufacturing. Berlin Heidelberg: Springer-Verlag; 2002. [5] Gursoy HN, Patrikalakis NM. Automated interrogation and adaptive subdivision of shape using medial axis transform. Advances in Engineering Software and Workstations 1991;13(5/6):287–302. [6] Elber G, Cohen E, Drake S. MATHSM: Medial axis transform toward high speed machining of pockets. Computer-Aided Design 2004;37:241–50. [7] Cao L, Liu J. Computation of medial axis and offset curves of curved boundaries in planar domain. Computer-Aided Design 2008;40:465–75. [8] Fabri R, Estrozi LF, Costa LF. On voronoi diagrams and medial axes. Journal of Mathematical Imaging and Vision 2002;17:27–40. [9] Žalik B. An efficient sweep-line Delaunay triangulation algorithm. ComputerAided Design 2005;37(10):1027–38. [10] Held M. Voronoi Diagrams and offset curves for curvilinear polygons. Computer-Aided Design 1998;30(4):287–300. [11] Chou JJ. Voronoi diagrams for planar shapes. IEEE Computer Graphics and Applications 1995;(March):52–9. [12] Farouki T, Ramamurthy R. Degenerate point/curve and curve/curve bisectors arising in medial axis computations for planar domains with curved boundaries. Computer-Aided Geometric Design 1998;15:615–35. [13] Ramamurthy R, Farouki T. Voronoi diagram and medial axis algorithm for planar domains with curved boundaries I: Theoretical foundations. Journal of Computational and Applied Mathematics 1999;102:119–41. [14] Ramamurthy R, Farouki T. Voronoi diagram and medial axis algorithm for planar domains with curved boundaries II: Detailed algorithm description. Journal of Computational and Applied Mathematics 1999;102:253–77. [15] Seong J.K., Cohen E., Elber G.. Voronoi diagram computations for planar NURBS curves, in: Symposium on Solid and Physical Modeling 2008; 67–77. [16] Pottmann H, Wallner J. Computational line geometry. Springer Verlag; 2001.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059 [17] Choi HI, Choi SW, Moon HP. Mathematical theory of medial axis transform. Pacific Journal of Mathematics 1997;181(1):57–88. [18] Choi HI, Choi SW, Moon HP, Wee NS. New algorithm for medial axis transform of plane domain. Graphical Models and Image Processing 1997;59(6):463–83. [19] Degen WLF. Exploiting curvatures to compute medial axis for domains with smooth boundary. Computer Aided Geometric Design 2004;21: 641–660. [20] Siddiqi K, Bouix S, Tannenbaum A, Zucker SW. Hamilton–Jacobi skeletons. International Journal of Computer Vision 2002;48(3):215–31. [21] Damon J. Geometry and medial structure. In: Pizer S, Siddiqi K, editors. Medial representations: Mathematics, algorithms, and applications. Springer/Kluwer series Comp. Imag. and Vision, vol. 37. Springer-Verlag; 2008.
1059
[22] Oliveira JB, De Figueiredo LH. Robust approximation of offsets, bisectors, and medial axes of plane curves. Reliable Computing 2003;9(2):161–75. [23] Chiang ChSh, Hoffmann ChM, Lynch RE. How to compute offset without selfintersections, In: Silbermann MJ, Tagare HD, editors. Curves and Surfaces in Computer Vision and Graphics II 1992;1610:76–87. [24] Choi HI, Choi SW, Han ChY, Kim T, Kwon SH, Moon HP, Roh KH, Wee NS. Twodimensional offsets and medial axis transform. Advances in Computational Math 2008;28:171–99. [25] Farin G. Curves and surfaces for computer aided geometric design. 5th ed. Morgan Kaufmann; 2002. [26] Berrut JP, Trefethen LlN. Barycentric Lagrange interpolation. SIAM Review 2004;46(3):501–17.
Contents lists available at ScienceDirect
Computer-Aided Design journal homepage: www.elsevier.com/locate/cad
Medial axis of a planar region by offset self-intersections R. Dorado ∗ Department of Mechanical and Mining Engineering, University of Jaén, Edificio A-3, Campus Las Lagunillas, 23071 Jaén, Spain
article
info
Article history: Received 31 December 2008 Accepted 22 August 2009 Keywords: Medial axis Medial axis transform Offset self-intersection Interpolation
abstract The medial axis (MA) of a planar region is the locus of those maximal disks contained within its boundary. This entity has many CAD/CAM applications. Approximations based on the Voronoi diagram are efficient for linear-arc boundaries, but such constructions are more difficult if the boundary is free. This paper proposes an algorithm for free-form boundaries that uses the relation between MA and offsets. It takes the curvature information from the boundary in order to find the self-intersections of successive offset curves. These self-intersection points belong to the MA and can be interpolated to obtain an approximation in Bézier form. This method also approximates the medial axis transform by using the offset distance to each self-intersection. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction
a The Medial Axis (MA) or skeleton of a planar region is the locus of centers of maximal disks contained within its boundary. Together with the function that gives the radius of each disk, the MA defines a Medial Axis Transform (MAT). This concept, from the differential geometry, was proposed to describe biological shapes by Blum [1, 2], who also described two definition models, a maximal disk model looking for a circle that touches two different boundary points, and a grass–fire model taking the region’s boundary as initial fire and propagates within the region until it intersects itself. The MAT properties are very well-known because of the numerous studies about this concept, mainly for planar regions. Ramanathan and Gurumoorthy [3] enumerated these properties. The MAT simplifies regional representation, and retains the original information as Patrikalakis and Maekawa [4] recalled. This is interesting in many computer animation applications because changing MA or the radius function is enough to deform the shape of a region, to construct finite element meshes from MA [5], or during numerical control machining when defining an optimum tool path for operations such as pocketing [6]. Despite the numerous efficient methods to construct the MA of domains with linear-arc boundaries, the MA of a free-form region is still a problem [7]. One option is to approximate the boundary using straight lines, and compute the MA by one of the existing methods for polygons. The drawback is the topological error of the solution because MA is sensitive to boundary continuity (Fig. 1). This error is a problem in applications where the planar region has to be reconstructed from its MA, e.g. to deform a region via the MA. Studies of related topics, such as the Voronoi diagram and the bisector, are other ways of computing this problem. The
∗
Corresponding author. Tel.: +34 953 212439; fax: +34 953 212870. E-mail address: [email protected].
0010-4485/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cad.2009.08.005
b
Fig. 1. MA computed from: (a) exact boundary (b) polygonal approximation.
MA of a free-boundary region can be computed using the Voronoi diagram [4,8], a dual of Delaunay triangulation [9]. Some scholars have studied the bisector to compute Voronoi diagrams. Held [10] for curvilinear polygons and Chou [11], Farouki and Ramamurthy [12], Ramamurthy and Farouki [13,14], and Seong, Cohen and Elber [15] for curved boundaries have obtained the Voronoi diagram by trimming bisectors. The previous methods use distance functions, as Seong et al. [15] that compute the bisector curve as an intersection of two distance surfaces, it is similar to the relation between untrimmed MA segments and the intersection between developable surfaces proposed by Pottman and Waller [16]. On the other hand, the trimming operation increases the process complexity. Moreover, if the boundary is free, neither the Voronoi diagram nor the MA is, in general, a subset of the other [13]. To obtain the MA you have to trim and/or add segments of the Voronoi diagram. Recently, scholars have been searching for other solutions. Choi, Choi and Moon [17] proposed a decomposition lemma to divide a general region into fundamental domains, where it is easy to compute the MA. They presented a computation method [18] based on the following idea; find the MA vertices, divide the domain, and compute normal MA points between two vertices with a method analogous to the maximal disk model.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
1051
b
a
Fig. 2. (a) MAT definition. (b) Criterion of minimum distance.
Ramathan and Gurumoorthy [3] used a tracing algorithm based on the evaluation of both distance and curvature criteria. Another tracing algorithm was proposed by Degen [19], who looked for the centers of maximal disks inside a planar region, and used the curvature information and the theory of dupin cyclides for his predictor/corrector algorithm. Based on the moving Frenet frames and Cesareo’s approach, Cao and Liu [7] discovered an interesting relation between two adjoining curves and the MA, and created a suitable tracing algorithm for solving the problem for free boundaries. Existing methods to compute the MA of a free-form region have, in general, a relatively low accuracy [7]. Moreover, most of them use a polygon approximation of the boundary, so the correction of topological errors is difficult to guarantee. An additional processing is also required to eliminate the artificial segments that may arise due to the discretization [3]. In this paper, the author proposes a method to approximate the MA of a free-boundary planar region. This method uses the exact analytic boundary, thus is free of the topologically incorrect segments produced by polygonal or curvilinear approximations. The solutions are MA segments and therefore, additional trimming operations, as in Voronoi techniques, are not needed. On the other hand, the approximation accuracy depends on the polynomial degree used to interpolate the computed MA points. To check the accuracy, a measure of the interpolation error is defined. The algorithm is analogous to the grass–fire model, because it finds the self-intersections of the successive boundary offset belonging to the MA. Note that grass–fire methods, like the proposed by Siddiqi et al. [20] or the methods reviewed by Damon [21], propagate a front through all pixels within a region of a digital image in order to choose possible candidates of MA, Oliveira and De Figueredo [22] developed a similar idea via a function distance. Consequently, it is difficult to guarantee that the solution is a thin set, i.e. without an interior, or is a connected set, and to compute the MA points accurately. By contrast, the proposed method uses a continuous definition of the boundary and its successive offsets instead of a digital image of the region, hence it guarantees an accurate, thin and connected (by interpolation) solution. Self-Intersection points are offset points with two or more different parameter values. For a parameter value u+ , it is possible to find a different parameter u− and the offset distance d for these parameter values, by solving a non-linear equation system. The solution for this system can be obtained using the Newton–Rapshon method with initial values of d and u− . These initial values can be approximated using information on the boundary curvature. In order to improve the algorithm and to deal with the least possible computations, this method interpolates the solutions at Chebyshev points. The solution is simple, natural, and thanks to the polynomial interpolation, with the same boundary
Fig. 3. Curvature sign and radius of curvature.
parameter. Finally, the algorithm gives the MAT, by interpolation of the offset self-intersection distance. This approximation is defined for regions with a general piecewise boundary of at least G1 (equal tangent direction at the piecewise unions). The extension to C 0 boundaries (equal points at the piecewise unions) is not the goal of this method, but Section 5.3 describes a simple example of this extension. The paper is arranged as follows. Section 2 describes the MA and offset relations. Section 3 shows a MA point classification and a definition of fundamental domains. It is possible to find the MA of general boundaries by the decomposition lemma [17]. If you can solve the fundamental case, you can solve the general one. Sections 4 and 5 describe the procedure for computing the MA points by offset self-intersections. The extension to general boundaries is based on the previous method and the minimum distance criterion, as shown in Section 6. Finally, some general examples and main conclusions are drawn in Sections 7 and 8 respectively. 2. Medial axis and offset self-intersections The MA or skeleton cs (u) of a planar region Ω can be expressed as a function of its boundary c(u), a radius function r (u) that gives the radius of maximal disks, and the unit normal directions from the boundary n(u). Let c(u), u ∈ [u0 , un ] be a parameterized closed boundary spline curve without intersections, in such a way that the region Ω is always to the left of parameter’s increasing direction, then the MA is given by (Fig. 2a): cs (u) = c(u) + r (u)n(u).
(1)
The unit normal n(u) and the curvature k(u) of c(u) are given by [4]:
(cu × cuu ) · ez , (2) |cu |3 where t(u) = cu / |cu | is the unit tangent of c(u), and ez = {0, 0, 1} n(u) = ez × t(u),
k(u) =
is a unit normal for the region plane. The positive curvature corresponds to c(u) points where the normal vector points out toward inside of the domain (Fig. 3).
1052
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
On the other hand, the MAT is the ordered set of all pair MA points and radius functions {cs (u), r (u)}. Note that MA and MAT depend on the radius function r (u). The Eq. (1) defines the centers of maximal disks. A disk is maximal when it is inside a region and does not contain other disks. According to Fig. 2b, the minimum distance dmin from the boundary to a disk center inside it and the disk radius R determine if the region contains the disk. These situations are possible and defining if a disk center belongs to the MA:
• dmin = R. Disks are within the region and have contact components (points or arcs) with the boundary. There are two cases: (i) One contact component. If the disk is an osculating one with the same tangent and curvature that the boundary (Fig. 3), its center is a MA point, because it is the maximal disk within the region with that contact (disk D1 in Fig. 2b). Otherwise, the disk is not maximal, so its center is not a MA point (disk D2 in Fig. 2b). (ii) More than one contact component. The centers are MA points, because they are the maximal disks within the region with those contact components (disk D3 in Fig. 2b). • dmin > R. This is the condition for all disks within the region without contacts with the boundary. They are not maximal, so their centers are not MA points (disk D4 in Fig. 2.b). • dmin < R. Disks intersect the boundary. They are not contained in the region and thus their centers are not MA points (disk D5 in Fig. 2b). The radius function r (u) defines the radii of maximal disks, and it is related to the offset distances. The proposed method uses this relation to compute the MA. In view of the applications for offset curves as tool path generation, feature recognition, or definition of tolerance regions, the previous idea is known, but it is used in the inverse problem to compute trimmed offsets. Chiang, Hoffmann and Lynch [23] use the MA to trim offsets, their idea is based on the fact that the MA cuts the offset of a closed curve at the self-intersections or at irregular offset points (isolated points and cusps). Choi et al. [24] show a completed study about the MA and offset relation. They note the coincidence between offset self-intersections and MA points, and use this relation to obtain a domain decomposition for offsets. Finally, Cao and Liu [7] also use the same idea to compute trimmed offset curves. The following definitions are useful to understand the relation between offset and MA. The offset cd (u) of a closed curve c(u) (with respect to distance d from the boundary to the interior of the region) is defined by: cd (u) = c(u) + dn(u).
(3)
The unit offset tangent td (u) and the offset curvature kd (u) depend on the tangent t(u) and the radius of curvature (radius of the osculating circle, Fig. 3) ρ(u) = 1/ |k(u)| of c(u): td ( u) = kd (u) =
1 − (d/ρ(u)) t(u), |1 − (d/ρ(u))|
(4)
1/ρ(u) . |1 − (d/ρ(u))|
(5)
Fig. 4. Relation between offset self-intersections and MA points.
Fig. 5. Decomposition lemma.
words, a self-intersection point is at the same distance respect to various boundary points and thus it belongs to the MA. Recall that the previous idea is used to trim offset curves in practical applications. This paper proposes the inverse problem, by computing the MA from offset self-intersections at different distances, as Fig. 4 shows. Instead of propagating the offset into the region, the algorithm solves a non-linear system of 2 equations for increasing parameter values u+ i , i = 0, 1, . . . , n. The solutions are the offset distances di to the MA and the parameter values u− i where the MA disks have a second contact point with the boundary. 3. Fundamental domains of a region Choi et al. [17,18] show that for a domain Ω with finite boundary curves linked with at least C 0 continuity (a spline curve), the MA is a connected geometric graph with a finite number of vertices and edges. Their decomposition lemma proves that a domain Ω can be divided into subdomains Ωj ,j = 0, 1, . . . , m, whose MAs are the Ω MA edges and vertices. A classification of MA points helps to obtain domain decomposition. According to Fig. 5, the MA points can be:
• Normal points. Center of maximal disks, of which each disk has two contact points with the boundary.
• Terminal points. These points are origin or end of one MA edge.
When distance d is equal to the radius of curvature, the offset has singularities (4) and (5). If d is larger than the minimum radius of curvature, the offset curve shows self-intersections. The Fig. 4 shows an example, note that ordinary MA points are centers of maximal disks with two contact points. The symbols +, − are used to distinguish between the parametric values of each point. Given the medial axis Eq. (1) and the offset Eq. (3), a self-intersection point with parameter values u+ , u− and offset distance d = r (u+ ) = r (u− ), lies on the medial axis. In other
(i) End-points are the centers of disks that have 1 contact component with the boundary: 1 point or 1 common circular arc. (ii) Corner points, boundary points where the unit tangent vector field is discontinuous. • Branch points. These points are the origin of more than 1 MA edges. (i) Non-generic points have 2 contact components with the boundary. (ii) Bifurcation points have 3 or more contact components with the boundary.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
1053
+ + boundary parameters u+ i , i = 0, 1, . . . , n, where u0 and un are the parameter values of the end-points, the corresponding di and u− i values to each MA point are the solution to the next selfintersection equations (+ and − are used to distinguish between the two contact components of the MA points): − cdi (u+ i ) = cdi (ui ),
(6)
the unknowns are di , ui . The non-linear system (6) can be solved by the Newton method, which has a quadratic convergence near a root (computational accuracy used < 10−8 ). The solutions u− i and + di for each u+ i give the MA points cs (ui ) by MA equation (1). Note − that to recover the boundary from the MA, u+ i and ui are necessary. The end-point parameters and radii of curvature give the initial values for solving the Eq. (6). The offset tangent and curvature equations (4) and (5) show that self-intersections begin when the offset distance is equal to the minimum radius of curvature ρ(u), thus self-intersection points start at the end-points computed in − Section 4.1. The MA computation begins with u+ 0 , u0 (for end+ + − points with 1 contact point u = u ), d0 = ρ(u0 ) = r (u+ 0 ) corresponding to one end-point. For u+ and the initial values i + − u− i−1 , di−1 = ρ(ui−1 ), the solutions for Eq. (6) are ui , di . Repeat this process until the parameter value corresponds to the other end+ + point u+ n , dn = ρ(un ) = r (un ). Finally, polynomial interpolation of MA points cs (u+ i ) and offset distances di = r (u+ ) gives the approximations of MA i and the radius function respectively in Bézier form, the standard use in practical applications [25]. This representation is easy to incorporate into CAD software, is easier to implement than the rational one, and admits a rational representation with unit weights. The main problem for obtaining an efficient interpolation is how to choose the u+ i values. A polynomial interpolation of equally spaced points (u+ i ) is highly ill-conditioned: small changes in the data may cause huge changes in the interpolant [26]. A numerically stable approach is to use point-sets that are clustered at the end-points of the interval with an asymptotic density −
Fig. 6. End-points computation.
In the following sections, the term vertices is used instead of terminal or/and branch points when the context refers to both types of MA point. A fundamental domain is a subdomain whose skeleton only has 2 vertices (terminal or/and branch points), and normal points. They are the smallest subdomains of a region and, thus, suitable cases for showing the performance of a computation method for normal MA points. 4. MAT computation for fundamental domains This section describes how to get the terminal points of a fundamental domain and how to use them to get its MA edge. For a C 0 boundary the terminal points are the boundary convex corners. Otherwise, a fundamental domain has end-points. 4.1. End-points computation Fundamental domains have no branch points. If the boundary is a parametric spline curve of at least G1 , the MA is a curve with two end-points. These points correspond to the disks with radii of curvature minimum within the region and therefore as different research works note [11,18], the boundary curvature allows for computing them. From the normal and curvature definition (2), and the minimum distance criterion, the end-points are the center of the osculating circles where the positive curvature k(u) attains its maxima. To obtain the end-points: 1. Compute the parameters’ values u where the curvature attains its maxima. If the curvature has a discontinuity at u, choose the u value of the greater curvature. Fig. 6 shows the parameters selected for 2 different end-points. 2. Compute the minimum curvature radii ρ(u) as the inverse of maximal curvatures. 3. Evaluate whether the disks with these minimum radii are within the domain. The minimum distance from the boundary to the center of the disk must be equal to its radius. As Cao and Liu [7] propose, the minimum distance parameter belongs to a set of parameters where the tangent cu (u) is perpendicular to the direction between the disk center p and the boundary c(u). The minimum distance associated with the roots of (p − c(u)) · cu (u) = 0 is the minimum distance needed. 4. The selected centers of the disks within the domain are the endpoints. Note that non-generic points are also MA points with maximum curvature. The number of contact components indicates if a point computed by the previous algorithm is non-generic. 4.2. MAT computation by offset self-intersections A self-intersection point is also a point of the MA, since the offset distance d, the parameter value u at this point, and the radius function are related r (u) = d. For a set of ordered
−1/2
proportional to 1 − u2 as n + 1 → ∞. The families of Chebyshev points are examples of these clustered point sets. These points are obtained by projecting equally-spaced points on the unit circle down to the unit interval [−1, 1]. Therefore, to compute a well-conditioned interpolation in a general domain u ∈ [u0 , un ], the u+ i can be the set of points: +
ui = 1 − cos
iπ n
+ u+ n − u0
2
+ u+ 0 ,
i = 0, 1, . . . , n.
(7)
The previous equation is a linear transformation of the Chebyshev points. Therefore, the point set (7) has the same density, between the parameters of two end-points, that the Chebyshev points over the unit interval. The exact MAT solution (MA and the radius function) is, in general, unknown, so the interpolation accuracyis also unknown. On the other hand, a valid interpolation MAT set c∗s (u), r ∗ (u) has
to verify Eq. (1): n(u) · c∗s (u) − c(u) /r ∗ (u) = 1, so a way of measuring the interpolation error s is:
n(u) · c∗s (u) − c(u) s = 1 − . r ∗ (u)
(8)
If s is larger than a desired value, we need more interpolation data points. A similar procedure is used to recover the boundary by substituting the MAT interpolation at Eq. (1). The boundary approximation error c (u) can be measured by the distance between the
1054
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
b
a
Fig. 7. MAT of an ellipse region (a) Degree 2 approximation (b) Degree 5 approximation.
Fig. 8. MAT approximation of a fundamental domain bounded by a degree-3 B-spline curve. (a) Degree 4 approximation (b) Degree 6 approximation.
boundary c(u) and the approximation c∗ (u). These curves have a common parameterization because of the approximation method.
c (u) = c(u) − c∗ (u) .
(9)
Steps to compute the MAT by offset self-intersections. 1. Compute the boundary end-points by the curvature information. For a general end-point with one contact component d0 = − + − ρ(u+ 0 ) = ρ(u0 ), and dn = ρ(un ) = ρ(un ). + 2. Choose the parameter values ui , i = 0, . . . , n with u+ 0 the parameter of one end-point and u+ n the other end-point parameter. Select a Chebyshev sequence (7).
3. Compute self-intersection points for each u+ i solving the Eq. (6). To compute the first self-intersection parameter u− 1 and d1 , + − − + + take u+ 1 and the initial values d1 = ρ(u0 ), u1 = u0 −(u1 − u0 ). The initial values for each new step are the solution of the + previous one until un . 4. Interpolate u+ , d i data points to obtain a Bézier MA and MAT i approximations. Finally, recover the boundary with the MAT set. In the following examples, the accuracy used to compute the curvature information and the self-intersections is 10−8 . The Figs. 7–14 were executed on a Mac with an Intel Core 2 Duo 2.8 GHz processor and 4GB RAM.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059 Table 1 CPU time used to compute Figs. 7–11.
a
1055
b
CPU time (s) Figure
End-points
MA edge
Total
7a 7b 8a 8b 9 10 11a 11b
0.018
0.008 0.020 0.025 0.061 0.010 0.030 0.010 0.070
0.026 0.038 0.055 0.091 0.220 0.290 0.030 0.270
0.030 0.210 0.260 0.020 0.200
Fig. 11. Degree 5 MAT approximation between different terminal points.
Fig. 9. Degree 2 MAT approximation between end-points with 1 contact arc. The boundary is a G1 spline curve.
End-points with 1 contact point have one parameter value. Therefore, the previous algorithm takes as initial values u+ 0 = u0 and u− = u . 0 0 A simple example is the MA of an ellipse, as shown in Fig. 7. Note that this method recaptures the exact solution for MA. On the other hand, the MAT and the recovered boundary are approximations. The previous result improves with more interpolation points. A general case is drawn at Fig. 8. The exact solution is unknown, so the interpolation accuracy can be measured by the interpolation error (8). The computational time is less than 0.1 s and depends on the interpolation degree (Table 1). The CPU time is mostly spent to compute the end-points. On the other hand, the number of interpolation points (degree interpolation + 1) increases the time to approximate the MA. 5.2. End-points with one contact circular arc A circular arc is common to those disks with center at these endpoints. The end-points are the MA’s for each circular arc. The initial − parameter values u+ 0 and u0 for computing the MA edge, are the arc terminal points. Fig. 9 is a simple example that shows how the offset collapses in the MA. There is a parallel pair of boundary segments in Fig. 9, note that Newton method gives one solution for each parameter value of the Chebyshev sequence along the boundary. Although the results are, in general, approximations, in this example the method provides an exact solution for MAT. The algorithm spends 0.22 s (Table 1) to compute a solution free of topological errors. 5.3. Domain with corners If a boundary region has convex corner points (with curvature
≥ 0 around them), the MA begins and/or ends at these points.
Fig. 10. Degree 5 MAT approximation for a boundary with convex corner points.
5. Examples of fundamental domains
For fundamental domains with 2 convex corner points, the initial − values for computing the MA are u+ 0 = u0 and u0 = u0 , where u0 is the parameter value of one corner. Fig. 10 shows an example. Note that concave corners are never on the MA, once they have been detected and isolated, the proposed algorithm can also be used in this case. Other examples of MA between different types of terminal points are drawn in Fig. 11. The CPU time is less than 0.3 s, and as the previous cases, the end-points consume almost all this time.
5.1. End-points with one contact point 6. Extension to general domains The previous method gives an approximation for different planar domains without branch points. The next examples show practical results and the extension to boundaries with corner points. Table 1 indicates the computational time.
Once the vertices (terminal or/and branch points) are known, the MA approximation of a general domain is the union of all fundamental MA computed between 2 successive vertices. The main
1056
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
Fig. 12. Degree 5 MAT approximation of a general region with 1 bifurcation point.
problem with general domains is computing the branch points, because the algorithm for fundamental domains helps to compute the end-points. Consider that we have two vertices and between them there is 1 unknown branch-point. Using the self-intersection algorithm between the first 2 vertices, some of the solutions are not MA points, because the minimum distance from the boundary to these points is less than the associated offset distance. The last correct MA point approximates the desired branch point. The previous idea helps to build an iteration procedure based on: firstly, compute an approximation between the 2 vertices and use the last correct MA point and the first wrong MA point to compute the next step, repeat the process until the difference between the 2 points computed in the last step is less than the desired precision, then take the last correct MA point as the branch
point. The results of this procedure are branch point coordinates and 2 parameter values where the associated disk has contact points with the boundary. The rest of the branch-point parameter values are where the boundary has the minimum distance with the branch-point. These steps summarize a self-intersecting algorithm for general domains that follows an analogous strategy to the Cao and Liu [7] one. 1. Define an ordered vector v with the couples of end and nongeneric parameters computed with the boundary curvature, as explained in Section 4.1. Sort the couples of parameters into smaller to greater values. The branch-point parameters are divided and ordered in couples from the first contact point until the same point again.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
1057
Table 2 CPU time used to compute the examples of Section 7. Figure
s <
CPU time (s) Vertices
12 13 14a
14b
14c
Fig. 13. Degree 20 MAT approximation of a general region with 2 bifurcation points.
2. Compute the self-intersections between the first 2 elements of v. Apply the minimum distance criterion. If their points do not belong to MA go to step 3. If all points belong to MA then go to step 4. 3. Compute the branch-point between the 2 first elements of v. Add the branch points’ parameters to v and compute the MA points between the 2 new first points of v. 4. Delete the 2 first elements of v and go to step 2. 5. Repeat from steps 2 to 4 until vector v is empty. 7. Examples of general domains The next examples help to understand the performance of the previous general algorithm. Fig. 12 shows the case of a domain with 3 end-points and 1 bifurcation point. To compute the MAT of this example:
• Generate a vector v with the end-points parameter values. = u− There are 3 end-points at u+ i , i = 0, 1, 2, so v = i + + + u0 , u1 , u2 . • Compute self-intersections between the 2 first elements of v. Compute the bifurcation point between the first 2 endpoints, as explained in Section 6. The bifurcation point has 3 contact components with the boundary ub0 , ub1 , ub2 . Join and sort the parameters to obtain a
+ + b b b new vector v = u+ 0 , u0 , u1 , u1 , u2 , u2 . In this example all MA segments can be computed from end-points with one contact, but if the construction begins in a bifurcation point, you need 2 initial parameters. • Compute the MA points by the offset self-intersections between the elements v and, finally, remove elements 1, 2 to obtain 1, 2 of + b b . a new v = u+ , u , u , u 1 1 2 2 There is no bifurcation point between the elements 1, 2 of v. • Compute the MA points between the elements 1, 2 of v. Next, b delete these elements to obtain a new v = u+ 2 , u2 . There is no bifurcation point between the elements 1, 2 of v, so repeat this step. • Compute the MA points between 1, 2 elements of v. Next, eliminate these elements to obtain a new v = {}. The vector v is empty.
A 2-branch point example obtained following the same general algorithm is drawn at Fig. 13. Finally, more difficult examples are drawn at Fig. 14. The procedure is the same for regions with holes inside (Fig. 14b), the main different is that the minimum distance criterion and the algorithm to compute the end-points also have to consider the holes. Fig. 14c shows the MA of a leaf computed from a photograph. Table 2 shows the CPU time used to compute the Figs. 12– 14. Note that branch and terminal points can be used to approximate the MA edges with different interpolation degrees and consequently their computational time is fixed. However, the
−3
11 × 10 8 × 10−7 4 × 10−2 6 × 10−4 2 × 10−7 10−2 2 × 10−4 5 × 10−7 10−2 10−4 3 × 10−7
1.00 1.90 3.71
4.41
20.80
MA edges
Total
0.05 0.81 0.17 0.41 2.52 0.22 0.51 3.83 0.84 1.99 14.61
1.05 2.71 3.88 4.12 6.23 4.63 4.92 8.24 21.64 22.79 35.41
number of interpolation points increases the CPU time, and the accuracy. The results are obtained in seconds for complex boundaries as in Fig. 14. The recent algorithm presented by Cao and Liu [7] spends minutes to compute the MA of regions with a similar number of end-branch points and MA edges. 8. Conclusions The method proposed approximates the MA of a free-form region using the boundary information. The algorithm has a general purpose. It can be of interest in tool-path optimization, and deformation of planar regions, because it is easy to program, easy to incorporate into CAD/CAM software, and it spends seconds to compute the examples shown in Sections 5 and 7. The MA points are computed by offset self-intersections. The relation between offset and MA is known, but it is, in general, used to trim offset curves [23,24] instead of computing the MA as in this work. Once the MA vertices (terminal or/and branch points) are known by the boundary curvature and the minimum distance criterion, their parameter and radius function provide the initial information needed to begin a tracing algorithm through the domain boundary parameters. This tracing algorithm is based on solving a non-linear system of two equations where the unknowns are the radius function and the second parameter associated with the computed MA point. Nonlinear systems are not a great drawback, because we know initial values near the solution. In order to obtain a more efficient approximation, the computed points are interpolated. The parameter values for the computed MA are the Chebyshev points, because they make the interpolation more stable. The solutions are curves in Bézier form being the standard in CAD/CAM, and have the same parameter as the boundary. The self-intersection method also gives the MAT of a region by interpolation of the radius function, and allows for recovering the boundary curve. The recovered boundary error is the distance in respect of the original curve, and the MA equation (1) allows for measuring the error of the MAT interpolation. The computational time depends on the boundary curve, because it determines the number of MA edges, vertices, and the interpolation degree needed to get accurate results. The CPU time used to compute the MA is less than 1 min in the examples (Table 2). Some differences between the proposed method and main existing constructions:
• Grass–fire methods are applied to digital images and thus it is difficult to guarantee thin and connected sets, and to compute accurate solutions. The proposed method provides an accurate and polynomial solution.
1058
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059
b
a
c
Fig. 14. MAT of general regions, the boundary is a degree-3 B-spline (a) Hand (b) Hand with a hole (c) Oak leaf.
• Voronoi methods need additional trimming operations because if the boundary is free the Voronoi diagram is, in general, not the MA. Offset self-intersections are directly the MA points. • Maximal disks methods uses the differential properties of the MA and the boundary. The accuracy depends on the initial values, the tracing step size and the iterative step size. A corrector or predictor algorithm [19,7] improve the accuracy. The proposed algorithm is based on the relation between the offset distance and the MAT radius function. The accuracy depends on the computation of vertices and MA points (computational accuracy < 10−8 ) and the interpolation degree. Additional algorithms to improve the accuracy are not needed. Acknowledgements This research is supported by the Spanish Ministerio de Innovación y Ciencia, under research grants JC2008-00107, DPI200609091, DPI2009-10078, co-financed by the ERDF (European Regional Development Fund). The author thanks Dr. Borut Žalik for several helpful discussions on medial axis. References [1] Blum H. A transformation for extracting new descriptors of shape. In: Models for perception of speech and visual form. Weinant Wathen-Dunn MIT Press; 1967.
[2] Blum H. Biological shape and visual science (Part I). Journal of Theoretical Biology 1973;38:205–87. [3] Ramanathan M, Gurumoorthy B. Constructing medial axis transform of planar domains with curved boundaries. Computer-Aided Design 2002;35: 619–632. [4] Patrikalakis NM, Maekawa T. Shape interrogation for computer aided design and manufacturing. Berlin Heidelberg: Springer-Verlag; 2002. [5] Gursoy HN, Patrikalakis NM. Automated interrogation and adaptive subdivision of shape using medial axis transform. Advances in Engineering Software and Workstations 1991;13(5/6):287–302. [6] Elber G, Cohen E, Drake S. MATHSM: Medial axis transform toward high speed machining of pockets. Computer-Aided Design 2004;37:241–50. [7] Cao L, Liu J. Computation of medial axis and offset curves of curved boundaries in planar domain. Computer-Aided Design 2008;40:465–75. [8] Fabri R, Estrozi LF, Costa LF. On voronoi diagrams and medial axes. Journal of Mathematical Imaging and Vision 2002;17:27–40. [9] Žalik B. An efficient sweep-line Delaunay triangulation algorithm. ComputerAided Design 2005;37(10):1027–38. [10] Held M. Voronoi Diagrams and offset curves for curvilinear polygons. Computer-Aided Design 1998;30(4):287–300. [11] Chou JJ. Voronoi diagrams for planar shapes. IEEE Computer Graphics and Applications 1995;(March):52–9. [12] Farouki T, Ramamurthy R. Degenerate point/curve and curve/curve bisectors arising in medial axis computations for planar domains with curved boundaries. Computer-Aided Geometric Design 1998;15:615–35. [13] Ramamurthy R, Farouki T. Voronoi diagram and medial axis algorithm for planar domains with curved boundaries I: Theoretical foundations. Journal of Computational and Applied Mathematics 1999;102:119–41. [14] Ramamurthy R, Farouki T. Voronoi diagram and medial axis algorithm for planar domains with curved boundaries II: Detailed algorithm description. Journal of Computational and Applied Mathematics 1999;102:253–77. [15] Seong J.K., Cohen E., Elber G.. Voronoi diagram computations for planar NURBS curves, in: Symposium on Solid and Physical Modeling 2008; 67–77. [16] Pottmann H, Wallner J. Computational line geometry. Springer Verlag; 2001.
R. Dorado / Computer-Aided Design 41 (2009) 1050–1059 [17] Choi HI, Choi SW, Moon HP. Mathematical theory of medial axis transform. Pacific Journal of Mathematics 1997;181(1):57–88. [18] Choi HI, Choi SW, Moon HP, Wee NS. New algorithm for medial axis transform of plane domain. Graphical Models and Image Processing 1997;59(6):463–83. [19] Degen WLF. Exploiting curvatures to compute medial axis for domains with smooth boundary. Computer Aided Geometric Design 2004;21: 641–660. [20] Siddiqi K, Bouix S, Tannenbaum A, Zucker SW. Hamilton–Jacobi skeletons. International Journal of Computer Vision 2002;48(3):215–31. [21] Damon J. Geometry and medial structure. In: Pizer S, Siddiqi K, editors. Medial representations: Mathematics, algorithms, and applications. Springer/Kluwer series Comp. Imag. and Vision, vol. 37. Springer-Verlag; 2008.
1059
[22] Oliveira JB, De Figueiredo LH. Robust approximation of offsets, bisectors, and medial axes of plane curves. Reliable Computing 2003;9(2):161–75. [23] Chiang ChSh, Hoffmann ChM, Lynch RE. How to compute offset without selfintersections, In: Silbermann MJ, Tagare HD, editors. Curves and Surfaces in Computer Vision and Graphics II 1992;1610:76–87. [24] Choi HI, Choi SW, Han ChY, Kim T, Kwon SH, Moon HP, Roh KH, Wee NS. Twodimensional offsets and medial axis transform. Advances in Computational Math 2008;28:171–99. [25] Farin G. Curves and surfaces for computer aided geometric design. 5th ed. Morgan Kaufmann; 2002. [26] Berrut JP, Trefethen LlN. Barycentric Lagrange interpolation. SIAM Review 2004;46(3):501–17.