Computation of medial axis and offset curves of curved boundaries in planar domain

Computer-Aided Design 40 (2008) 465–475 www.elsevier.com/locate/cad

Computation of medial axis and offset curves of curved boundaries in planar domain Lixin Cao ∗ , Jian Liu School of Mechanical Engineering, Dalian University of Technology, China Received 10 October 2007; accepted 8 January 2008

Abstract In this paper, we begin our research from the generating theory of the medial axis. The normal equidistant mapping relationships between two boundaries and its medial axis have been proposed based on the moving Frenet frames and Cesaro’s approach of the differential geometry. Two pairs of adjoint curves have been formed and the geometrical model of the medial axis transform of the planar domains with curved boundaries has been established. The relations of position mapping, scale transform and differential invariants between the curved boundaries and the medial axis have been investigated. Based on this model, a tracing algorithm for the computation of the medial axis has been generated. In order to get the accurate medial axis and branch points, a Two Tangent Points Circle algorithm and a Three Tangent Points Circle algorithm have been generated, which use the results of the tracing algorithm as the initial values to make the iterative process effective. These algorithms can be used for the computation of the medial axis effectively and accurately. Based on the medial axis transform and the envelope theory, the trimmed offset curves of curved boundaries have been investigated. Several numerical examples are given at the end of the paper. c 2008 Elsevier Ltd. All rights reserved.

Keywords: Medial axis transform; Cesaro’s approach; Tracing algorithm; Iterative algorithm; Offset curve

1. Introduction The medial axis transform (MAT), also called skeleton, was first introduced by Blum [1] as a description of shape. Its classical definitions include grass-fire model and maximal disk model. The grass-fire model means that an object’s boundary is taken as an initial fire front that propagates within the object’s interior region. Points where the fire front folds or interacts with itself are retained as the skeleton points. The maximal disk model means that the medial axis of a planar domain is the locus of the center of a maximal disc, which touches the boundary in at least two points. The MAT of a planar domain is the medial axis (MA) together with its associated radius function. Since its introduction, the MAT has been paid great attention by many researchers. This is because the MAT possesses the following excellent characteristics, such as uniqueness, invertibility, symmetry, topological equivalence, one-to-one correspondence, etc. In other words, there is a unique MAT for a boundary curve and the boundary curve ∗ Corresponding author. Tel.: +86 0411 84708411.

E-mail address: [email protected] (L. Cao). c 2008 Elsevier Ltd. All rights reserved. 0010-4485/$ - see front matter doi:10.1016/j.cad.2008.01.002

can be reconstructed by its MAT uniquely. Based on these characteristics and researchers’ efforts, the MAT has been used in pattern analysis, image analysis, finite element mesh generation, CNC tool-path planning, etc. And it can potentially be used as a representation scheme in geometric modelers, along with constructive solid geometry (CSG) and boundary representation (B-rep) [2]. The algorithms for the MAT have been investigated for many years. According to different geometric objects, present algorithms can be classified as discrete MAT and continuous MAT [3,4]. For the discrete MAT, the object is stored in binary image. The commonly used algorithms for discrete MAT include thinning method and distance transformation method based on the grass-fire model. The thinning method mostly needs iteration processes which is time consuming. Furthermore, the generated result of the thinning method is not accurate enough due to the discrete processing. For a two-dimensional discrete domain, the traveling direction of any point is along one of 8-connected pixels at least. It cannot simulate the arbitrary traveling direction of the grassfire model exactly. The distance transformation method can generate skeletons with accurate position, but connectivity

466

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

of the skeleton cannot be guaranteed due to the discrete processing. In addition, the skeleton point based on the distance transformation is checked independently, the neighboring points are not considered mostly. This also results in the medial axis losing the characteristics of connectivity and 1-pixel width. There is no effective approach for this problem up to now. This kind of discrete MAT methods are mainly used in pattern analysis and image analysis, but not suitable for CNC tool path planning of a region with continuous boundary curves. The continuous MAT has played extensive roles in solid modeling and computational geometry. However, the application of the continuous MAT is restricted in some degree as yet. The main reason is that the algorithms to compute the MAT of a planar domain appeared slow and not accurate [5]. Before 1989, these algorithms are restricted to very simple cases like polygons or domains whose boundaries are composed by line segments and circular arcs. For domains with curved boundary, its medial axis is nonrational and the exact computation for MAT is very difficult. One simple approach is approximating the curved boundary by straight line segments in a given tolerance firstly, and then invokes the mature algorithms of MAT for polygons. This will result in wrong topological results, and the final results are influenced by the given tolerance. The essential reason is that the MAT is sensitive to the continuous order of the boundary. As a result, for a MAT of a domain with curved boundary, one must adopt exact representation of curved boundary to avoid wrong topological result. During the decade after 1989, some researchers began to study the curve/point bisector. Meanwhile, the main discovery is that the bisector of two curves can be considered as the envelope of a one parameter family of curve/point bisectors where the point moves on the second curve [6,7]. An algorithm considering the degeneracy was developed by Farouki [8] based on these investigations, such as point/curve bisectors where the point lies on the curve, curve/curve bisectors where the two curves are identical, and curve/curve bisectors for distinct curves that share a common endpoint. During this period, another recommendable algorithm is introduced by Choi [9] based on the domain decomposition lemma, which enables one to break up the complicated domain into smaller and simpler pieces where no branching point is included; thereby complicated processes for branching point can be avoided. However, the implementation of the decomposition algorithm is not easy, it needs multiple iteration processes. This limits the application of the decomposition algorithm greatly. In recent years, some theoretical progresses were made by the work of the following authors. Teixeira [10], August, Tannebaum and Zucker [11] investigated the relations of curvature, tangent and normal vectors between the boundary curve and its medial axis, a differential equation that gives the movement of the medial axis of a curve as this curve evolves through mean curvature motion has been established. Ramanathan [12] generated a tracing algorithm of MAT. The tracing procedure originates at a convex vertex, and then locates the next footpoint on one edge by the tracing step size and finds the corresponding footpoint on another edge enforcing

the distance criterion, thereby obtains the next medial axis point followed by checking against curvature criterion. This step requires the intersection of normals to the two edges and needs iteration processes to find the roots of a nonlinear equation in one variable. Thus, this algorithm is not a strict tracing algorithm. Emphasized on the boundary’s curvatures, Degen [5] developed a predictor/corrector algorithm by means of the theory of Dupin cyclides and the accuracy of computation of MA can be improved considerably. From the recent references, we know that the medial axis for the planar domains with curved boundary has been paid to more and more attentions by many scholars. At the same time, effective algorithms of MAT about the curved boundary have become a bottleneck for the application of MAT in engineering. Existing algorithms of MAT mostly need polygon approximation or numerical iteration process, correct topological result is difficult to guarantee. The computational accuracy and efficiency are relatively lower. Furthermore, deeper geometric properties of MAT have been exploited rarely, essentially only tangents come into play. In this paper, we begin our research from the generating theory of the medial axis. The normal equidistant mapping relationships between the two boundaries and its medial axis have been proposed based on the moving Frenet frames and Cesaro’s approach of the differential geometry [13]. The two pairs of adjoint curves have been formed and the geometrical model of the medial axis transform of the planar domains with curved boundaries has been established. The relations of position mapping, scale transform and differential invariants between the curved boundaries and the medial axis have been investigated. Based on this model, a tracing algorithm for the computation of the medial axis has been generated. In order to get the accurate medial axis and branch points of the curved boundaries, a Two Tangent Points Circle algorithm and a Three Tangent Points Circle algorithm have been generated, which use the results of the tracing algorithm as the initial values to make the iterative process effective. These algorithms can be used for the computation of the medial axis effectively and accurately. Based on the medial axis transform and the envelope theory, the trimmed offset curves of curved boundaries have been investigated. The global and local selfintersection points of the offset curves have been trimmed. Several numerical examples are given at the end of the paper. 2. Concomitant relations between medial axis and boundary curves In this section, we will start our research from the generating theory of the medial axis. Two boundaries and its corresponding medial axis have been regarded as two pairs of adjoint curves. Therein, the relations of position mapping, scale transform and differential invariants between the curved boundaries and the medial axis have been investigated based on the moving Frenet frames and Cesaro’s approach of the differential geometry [13]. The deeper research of the three groups of relations will provide theoretical foundation for the following tracing algorithm,

467

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

points q1 and q2 , respectively; r denotes the radius of medial axis circle at point p. The relations of position mapping of the boundary curves and the medial axis have been given out by Eqs. (1)–(3). 2.2. Relations of scale transform In order to establish the parametric differential mapping relations between the boundary curves and the medial axis, Eq. (3) has been differentiated as follows ) (1) (1) dr1 = dr + e2 dr + r de2 . (4) (2) (2) dr2 = dr − e2 dr − r de2

Fig. 1. Medial axis and boundary curves.

which is used to compute the medial axis of the curved boundaries. Suppose C1 and C2 are two boundary curves in a given plane, the set of the centers of circles that are tangent to C1 and C2 respectively at two distinct points is defined as medial axis C. The point p is called normal point of MA due to the circle touches exactly two separate boundary curves. The circle is called medial axis circle with radius r . The tangent points are q1 and q2 , as shown in Fig. 1. Frame {p, e1 e2 } is moving Frenet frame of the MA with its origin at point p, and frames (1) (1) (2) (2) {q1 , e1 e2 } and {q2 , e1 e2 } are moving Frenet frames of the boundaries C1 and C2 , with their origins at points q1 and q2 , respectively. s, s1 and s2 are natural parameters of C, C1 and C2 , respectively. 2.1. Relations of position mapping As shown in Fig. 1, the mapping relations of the Frenet frames of the curved boundaries and the medial axis can be expressed as: ) (1) e1 = e1 sin θ − e2 cos θ (a) (1) (1) e2 = e1 cos θ + e2 sin θ (b) ) (2) e1 = e1 sin θ + e2 cos θ (a) (2) (2) e2 = −e1 cos θ + e2 sin θ (b) where θ denotes the angle between radius vector q1 p or q2 p and the unit vector e1 at point p. If we regard the two boundary curves as two adjoint curves of the medial axis, the boundary curves can be written as the following equations considering the symmetric property of medial axis. ) (1) r1 = r + r e2 (3) (2) r2 = r − r e2 where r denotes radius vector of the medial axis at point p; r1 , r2 denote radius vectors of the boundaries C1 and C2 at

By using the Frenet equations of the planar curve of the differential geometry [13], Eq. (4) can be simplified as follows ) (1) (1) e1 (1 + r k1 ) ds1 = e1 ds + e2 dr (a) (5) (2) (2) e1 (1 − r k2 ) ds2 = e1 ds − e2 dr (b) where k1 , k2 denote the curvatures of boundary curves C1 and C2 at the points q1 and q2 , respectively. Computing dot (1) products of both sides of Eq. (5)(a) and (b) by using e1 and (2) e1 respectively, we have  sin θ  ds1  =  ds 1 + r k1 . (6) sin θ  ds2   = ds 1 − r k2 The parametric differential mapping relations between the medial axis and the boundary curves have been given out by Eq. (6). In the same way, computing dot products of both sides (1) (2) of Eq. (5)(a) and (b) by using e2 and e2 respectively, we have cos θ ds + dr = 0.

(7)

Eq. (7) gives the differential relation that r and θ should abide. 2.3. Relations of curvatures Differentiating Eqs. (1)(a) and (2)(a), we have (1)

de1 = de1 sin θ + e1 cos θdθ − de2 cos θ + e2 sin θdθ (2)

de1 = de1 sin θ + e1 cos θdθ + de2 cos θ − e2 sin θdθ

) . (8)

By using the Frenet equations of the planar curve of the differential geometry [13] and considering Eqs. (1) and (2), Eq. (8) can be simplified as follows  k1 ds1 = kds + dθ . (9) k2 ds2 = kds − dθ Substituting Eq. (6) into Eq. (9), the curvature of the medial axis can be solved as k=

sin θ (k1 + k2 ) sin θ (ρ1 + ρ2 ) = 2 (1 + k1r ) (1 − k2r ) 2 (ρ1 + r ) (ρ2 − r )

(10)

where ρ1 = 1/k1 , ρ2 = 1/k2 denote the curvature radii of the boundary curves C1 and C2 at the points q1 and q2 , respectively.

468

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

From Eq. (10), we know that this equation is meaningful if 1 + r k1 6= 0 and 1 − r k2 6= 0. It means that the normal point of the medial axis will degenerate into an end point when the radius of the MA circle is exactly equal to the radius of the convex boundary’s curvature at the tangent points. Hence, the conditions 1 + r k1 6= 0 and 1 − r k2 6= 0 can be used as a terminating condition for the following tracing algorithm. 3. Tracing algorithm of the medial axis For the given boundary curves C1 and C2 , in order to construct the tracing algorithm of the MA, an initial point p of MA and the corresponding tangent points q1 and q2 should be given. It can be acquired

by using iterative algorithm [12]. Then e2 = q1 q2 / q1 q2 2 , e1 is perpendicular to e2 , e3 points at reader and abide right hand rule. M denotes the curvature circle of medial

axis at point p, as shown in Fig. 2.
And r = p − q1 2 , θ = arccos (x1 /r ) = arccos q1 p • e1 /r . By using Eq. (10), we can get the curvature of the MA at point p. The vector equation of the next point p∗ of MA can be written as follows based on the Taylor’s expansion (further terms of order ≥3 are omited) 1 r p∗ = r (s) + 1se1 + k1s 2 e2 (11) 2 where 1s denotes the infinitesimal arc-length along e1 at the point p of the MA. The Frenet frame corresponding to the point p∗ can be updated as follows  e∗1 = e1 cos δ + e2 sin δ (12) e∗2 = −e1 sin δ + e2 cos δ where δ = 1s/ρ denotes an angle corresponding to 1s of curvature circle M, as shown in Fig. 2; ρ = 1/k denotes the radius of the curvature circle M. Based on Eq. (6), the corresponding arc-length parametric increments 1s1 and 1s2 of the boundary curves C1 and C2 can be acquired. The parametric increments of the boundary

curves C1 and C2 can be expressed as 1u 1 = 1s1 / r01 2 and 1u 2 = 1s2 / r02 2 , respectively. If we let  u ∗1 = u 1 + 1u 1 (13) u ∗2 = u 2 + 1u 2 and we put the parameters u ∗1 and u ∗2 into equations of the boundary curves C1 and C2 , we can get the new tangent points q∗1 and q∗2 which corresponding to the point p∗ . And the corresponding curvatures of the boundary curves at points q∗1 and q∗2 can be acquired by using the following formula.

0

r × r00 i i 2 ∗ k i = 3 i = 1, 2. (14)

r0 i 2 Thus, we can get the next point p∗ of MA, corresponding Frenet frame and the new tangent points q∗1 and q∗2 . Repeat this tracing process, all the points of MA and the corresponding radii of the medial axis circles can be obtained. Here we will summarize this algorithm as a pseudo code.

Fig. 2. Tracing algorithm of the medial axis.

Procedure Tracing Algorithm (p, q1 , q2 , 1s) do { Compute

unit vector e1 and e2 ; r := p − q1 2 ;
θ := arccos q1 p • e1 /r ; Compute curvature of MA at point p by using Eq. (10); Compute next point p∗ of MA by using Eq. (11); Update Frenet frame of MA by using Eq. (12); Compute arc-length parametric increments 1s1 and 1s2 of the boundary by using Eq. (6);

0curves

r ; 1u 1 := 1s1 / 1 2 1u 2 := 1s2 / r02 2 ; Compute new tangent points q∗1 , q∗2 and boundaries’ curvatures by using Eqs. (13) and (14); p := p∗ ; q1 := q∗1 ; q2 := q∗2 ; }while {1 + r k1 6= 0 and 1 − r k2 6= 0}; End Tracing Algorithm

4. Accurate algorithms of medial axis 4.1. Two Tangent Points Circle algorithm The above tracing algorithm can be used to compute medial axis of the curved boundaries of planar domain, but the computational accuracy of the tracing algorithm depends on the arc-length parametric increments 1s of the medial axis. If 1s is too large, the computational accuracy cannot satisfy the needs of the engineering application. In this section, we will investigate a Two Tangent Points Circle algorithm, which means a circle is tangent to the boundary curves C1 and C2 at two distinct points correctly. In order to make this algorithm effective, we use the results of the tracing algorithm as the initial values of this Two Tangent Points Circle algorithm. As shown in Fig. 3, point p of MA and the corresponding tangent points q1 (u 1 ) and q2 (u 2 ) of the boundary curves C1 and C2 are given by using

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

469

between the medial axis and the boundary curves have been built, then the accurate solutions can be acquired via less iterative processes. 4.2. Three Tangent Points Circle algorithm

Fig. 3. Two Tangents Points Circle algorithm.

the tracing algorithm. e1 and e2 are unit tangent vector and unit normal vector of the medial axis at point p, respectively. The minimum distances (d1 and d2 ) between the point p and the boundary curves C1 and C2 can be acquired by using an algorithm (see Section 4.3), which use the points q1 and q2 as the initial values. The points b1 and b2 are the corresponding footpoints of the boundary curves C1 and C2 , respectively. If d1 is larger than d2 , then let the point p moves a step 1 along the direction of e2 ; otherwise, let the point p moves a step 1 along the opposite direction of e2 . By decreasing the step 1 at the right place, the minimum distances d1 and d2 will be equal eventually. That means the accurate medial axis circle with center p∗ is acquired, as shown in Fig. 3. Here we will summarize this algorithm as a pseudo code. Procedure Two Tangent Points Circle (p, u 1 , u 2 ) Compute unit vector e1 and e2 ; d1 = MinDisBetweenPointAndCurve(p, u 1 ); d2 = MinDisBetweenPointAndCurve(p, u 2 ); If (d1 < d2 ) {reverse the direction of unit vector e2 ; switch d1 and d2 } do { If (d1 < d2 ) do { 1 := 1/2; p := p − 1e2 ; d1 = MinDisBetweenPointAndCurve(p, u 1 ); d2 = MinDisBetweenPointAndCurve(p, u 2 ); }while (d1 < d2 ); else { p := p + 1e2 ; d1 = MinDisBetweenPointAndCurve(p, u 1 ); d2 = MinDisBetweenPointAndCurve(p, u 2 ); } }while (|d1 − d2 | >tolerance); End Two Tangent Points Circle Since this algorithm is based on the results of the tracing algorithm, the parametric differential mapping relations

For a connected domain in the Euclidean plane, the medial axis is consisting of branch points, end points and the normal points generally [3–5,7]. A point whose medial axis circle is tangent to the domain boundary in three or more distinct points is called branch point, as shown in Fig. 4. In this section, we will investigate a Three Tangent Points Circle algorithm, which means a circle is tangent to the boundary curves at three distinct points correctly. During the process of the tracing algorithm or Two Tangent Points Circle algorithm, we should calculate the minimum distance of the current point p of medial axis from the boundary curves and check if the minimum distance is less than the radius of the current medial axis circle. If the distance is larger than the radius of the current medial axis circle, which means that the branch point has not been appeared, then continue the tracing algorithm and Two Tangent Points Circle algorithm. On the contrary, if the minimum distance is smaller than the radius of the current medial axis circle, which means the current point p has already gone beyond a branch point. The branch point can be acquired by retracting the current point p. Based on the Two Tangent Points Circle algorithm, it is easy to give out a pseudo code of the Three Tangent Points Circle algorithm as follows. Procedure Three Tangent Points Circle (p, u 1 , u 2 , u 3 ) Compute unit vector e1 of medial axis at point p; d1 = MinDisBetweenPointAndCurve(p, u 1 ); d3 = MinDisBetweenPointAndCurve(p, u 3 ); do { If (d3 > d1 ) do { 1 := 1/2; p := p + 1e1 ; d1 = Two Tangent Points Circle (p, u 1 , u 2 ); d3 = MinDisBetweenPointAndCurve(p, u 3 ); }while (d3 > d1 ); else { p := p − 1e1 ; d1 = Two Tangent Points Circle (p, u 1 , u 2 ); d3 = MinDisBetweenPointAndCurve(p, u 3 ); } }while (|d1 − d3 | >tolerance); End Three Tangent Points Circle As shown in Fig. 4, the point p∗ is a branch point which is tangent to boundary curves at three distinct points q∗1 , q∗2 and q3 . Now, replacing the point q∗2 by the point q3 , and invoke the tracing algorithm and Two Tangent Points Circle algorithm to find medial axis C11 . Once the medial axis C11 is finished, then go back to the point p∗ , and use the points q∗2 and q3 as the tangent points, to find another medial axis C12 .

470

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

Fig. 5. Minimum distance between a point and.

5. Main algorithm

Fig. 4. Branch point.

4.3. Minimum distance between a point and a curve In Sections 4.2 and 4.3, we have mentioned the algorithm of the minimum distance between a point and a curve. The minimum distance between a point and curve/surface has been researched by many scholars. In order to make this paper more readable, the minimum distance algorithm will be rewritten in the following based on the references [14,15]. Suppose point p of MA and the corresponding tangent point q of the boundary curves C (no singular point is contained) are given by using the tracing algorithm. From Fig. 5 we know, if the curve does not contain a singularity point, then have ru 1u + nd = 1r

(15)

where ru denotes tangent vector of curve C at the point q; n denotes the unit normal vector of the curve C at the point q; 1r denotes the vector from point q to the point p; d denotes the minimum distance between the point p and the curve C. Computing dot products of both sides of Eq. (15) by using ru , we have 1u = (1r • ru ) / (ru • ru ) .

(16)

Let u = u + 1u, and repeat the former process until |1u| ≤ ε, (ε denotes the computational precision). Then we can get the minimum distance between point p and the curve C. It can be expressed as d = 1r • n. The pseudo code of this algorithm is given as follows Procedure MinDisBetweenPointAndCurve(p, u) do { r := r (u) ; ru := ru (u) ; 1r := p − r; d = |1r| ; 1u := (1r • ru ) / (ru • ru ) ; u := u + 1u; }while (|1u| > ε); End MinDisBetweenPointAndCurve

(17)

We now combine the procedures developed so far into a single coherent procedure. The first step of the main algorithm is to compute the end points of medial axis by seeking the maxima of boundary’s curvature where it has positive value. All these end points form an initial end points table. Afterwards, running the tracing and Two Tangent Points Circle algorithm from one of the initial end points until another end point is emerged. During this process, the minimum distance between the current point P of medial axis and the discrete boundary curves should be calculated. If the minimum distance is larger than the radius of the current medial axis circle, that means the branch point has not appeared, then continue the tracing and Two Tangent Points Circle algorithm. On the contrary, if the minimum distance is smaller than the radius of the current medial axis circle, which means the current point P has already gone beyond a branch point. In this case, we should switch to the Three Tangent Points Circle algorithm to get the branch point p∗ . Furthermore, adding points p∗ , q3 and q2 to the initial end points table and replacing the current tangent point q2 by new tangent point q3 in the same time. Then, continue the tracing and Two Tangent Points Circle algorithm. After each tracing process, the beginning and ending points should be canceled from the initial end points table. The main algorithm will be terminated when the initial end points table is empty. Here we summarize our main algorithm as a pseudo code. Begin Main Algorithm Compute initial end points table of MA; do{ Run Tracing Algorithm; Run Two Tangent Points Circle; Compute the minimum distance dmin of the point P from the discrete boundary curves; If (dmin > r ) continue; else { Run Three Tangent Points Circle for branch point p∗ and tangent points q1 , q2 and q3 ; Add points p∗ , q3 and q2 to the initial end points table; Use points p∗ , q3 and q1 as the current initial end points; Update unit vector e1 and e2 ; } while (the initial end points table is not empty); End Main Algorithm

471

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

Rϕ = r E1 (ϕ)

(21)

where E1 (ϕ) = − sin ϕi + cos ϕj. Substituting Eqs. (20) and (21) into the following envelope condition
k × Rϕ • Rt = 0. (22) We have x 0 cos ϕ + y 0 sin ϕ + r 0 = 0.

(23) q

By dividing the both sides of Eq. (23) by x 0 2 + y 0 2 . And let q q 2 2 0 0 0 0 cos α = x / x + y , sin α = y / x 0 2 + y 0 2 , we have cos (ϕ − α) = − q Fig. 6. Untrimmed offset curve.

Offset, as a geometric operation, is widely used in CAD/CAM. For example, offset curves can be used as the tool paths of CNC machines. For a given curve C : r (t) = (x (t) , y (t)) in a plane, its offset curves can be written as

= x (t) − d q

x 0 2 (t) +

 +  y (t) + d q



y 0 (t) y 0 2 (t)

i 

x 0 (t) x 0 2 (t) +

y 0 2 (t)

j

(18)

where n (t) denotes the unit normal vector of the curve C at the given point; t denotes the parameter of the curve C and d denotes the offset distance; i and j are unit vectors of axes X and Y , respectively. This formula gives the untrimmed offsets, as shown in Fig. 6. The untrimmed offset may have local and global selfintersection which should be excluded for the application, such as tool path of CNC machining. Generally, the trimming process is time-consuming. Once the MAT is gotten, the trimming process becomes much easier. Based on the invertibility of the MAT, if we have MAT, the original boundary curves and the corresponding offset curves can be reconstructed from the envelope of the medial axis circles. Let r (t) = (x (t) , y (t) , r (t)) be the MAT of a planar domain, the formula of the family of the medial axis circles, i.e. generating circles, can be written as R = r (t) + r e(ϕ)

(19)

where e (ϕ) = cos ϕi + sin ϕj and ϕ denotes the parameter of the generating circle. Differentiating Eq. (19) with respect to t and ϕ respectively, we have Rt = x 0 i + y 0 j + r 0 e (ϕ)

x 02 + y02

.

(24)

Solving Eq. (24) for ϕ, we have

6. Offset curves and medial axis transform

rd = r (t) + dn (t) 

r0

(20)

ϕ = α + π ± arccos q

r0 x 02 + y02

.

(25)

This formula is the detailed envelope condition that the envelope curves should abide. Substituting Eq. (25) into Eq. (19), we have the envelope formula of the medial axis circles q    −x 0r 0 ∓ y 0 x 0 2 + y 0 2 − r 0 2  i R = x + r  x 02 + y02 q    −y 0r 0 ± x 0 x 0 2 + y 0 2 − r 0 2  j. (26) + y + r  x 02 + y02 If we use the arc-length as the parameter of the medial axis and considering Eq. (7), Eq. (26) can be simplified as the follows R = r (s) + r e(α ± θ )

(27)

where e (α + θ) = cos (α + θ) i+sin (α + θ) j denotes the unit normal vector of the boundary curve at the footpoint; α denotes the angle between the axis i and the unit vector e1 at the point p;θ is defined in Section 2.1, as shown in Fig. 1. Eq. (27) implies that the footpoint, i.e. tangent point between the medial axis circle and the boundary curve, is just the point which forms the envelope curve. This point is called characteristic point in differential geometry [13]. It is possessed by both the medial axis circle and the envelope curve, and with common normal vector on it. The offset curves of the boundary curves can be expressed as Rd = r (s) + (r − d) e (α ± θ)

(28)

where d denotes the offset distance. The parts of the offset curves where the corresponding radius of the medial axis circle r is less than d have to be trimmed in order to avoid the selfintersection of the offset curves. This trimming process becomes simpler.

472

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

Fig. 7. Digital photo and boundary points.

Fig. 9. MA and medial axis circles.

Fig. 8. MA and branch points.

Fig. 10. MA and offset curves.

7. Numerical examples In the following, we will verify the above theoretical analysis and algorithms by means of four examples. The initial points of the boundary curves are acquired by using processes of digital photoing, LP filtering, HP filtering, edge enhancing, etc. The boundary curves of these examples are the cubic uniform Bspline curve. The control points of the B-spline curves were acquired by using the fair fitting method with the least squares of the B-spline curve. The initial points of the tracing algorithm can be found by seeking the maxima of boundary’s curvature where it has positive values. The computational accuracy or tolerance of the Two Tangent Points Circle algorithm, Three Tangent Points Circle algorithm and minimum distance algorithm is 10−7 . The computational time (6 min for the first example, 10 min for the second example, 12 min for the third example and 7 min for last example) of my examples is mainly

occupied by the computing process of the minimum distance of the point P of MA from the discrete boundary curves at each step in order to find the branch point. The extraction of the branch points of medial axis is a very difficult problem and has not been solved thoroughly up to now from the knowledge of the literatures. The results of these examples are given as follows. Example 1. Leaf of liriodendron. See Figs. 7–10. Example 2. Leaf of phoenix tree. See Figs. 11–14. Example 3. Leaf of oak. See Figs. 15–18. Example 4. Leaf of liriodendron with two holes. See Figs. 19– 21.

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

Fig. 11. Digital photo and boundary points.

473

Fig. 14. MA and offset curves.

Fig. 12. MA and branch points.

Fig. 15. Digital photo and boundary points.

The above three examples are all the star-shape topology. In order to verify the effectiveness of our algorithms, another example which has different topology from the above examples has been given in the following. 8. Conclusions

Fig. 13. MA and medial axis circles.

• By using the moving Frenet frames and Cesaro’s approach of the differential geometry, the normal equidistant mapping relationships between the two boundaries and its medial axis have been established. The two pairs of adjoint curves have been formed and the geometrical model of the medial axis

474

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

Fig. 16. MA and branch points.

Fig. 18. MA and offset curves.

Fig. 19. MA and branch points.

Fig. 17. MA and medial axis circles.

transform of the planar domains with curved boundaries has been generated. • The relations of position mapping, scale transform and the differential invariants between the curved boundaries and the medial axis have been investigated. Based on this model, a tracing algorithm for the computation of the medial axis has been generated. This tracing algorithm doesn’t need iteration, so it can be used for the computation of the medial axis effectively. By using this algorithm, the next point of the MA is acquired on the curvature circle

of the MA at the current MA point, the omitted errors are only third order infinitesimal. At the same time, the next footpoints are acquired by means of exact boundaries. This will have a correction for the results of the tracing algorithm. • Based on the tracing algorithm of the medial axis, two iterative algorithms have been generated. These algorithms can be used to compute medial axis accurately and find the branch points. This is very helpful for improving the computational accuracy of the tracing algorithm. The numerical examples also indicate that the accuracy of these algorithms is recommendable. By using the medial axis transform and the envelope theory, the trimmed offset curves of curved boundaries have been investigated.

L. Cao, J. Liu / Computer-Aided Design 40 (2008) 465–475

475

acknowledged. In addition, the authors would like to thank the anonymous reviewers for their considerate comments and suggestions. References

Fig. 20. MA and medial axis circles.

Fig. 21. MA and offset curves.

Acknowledgements This research is supported by the National Natural Science Foundation of China (50775022). This support is gratefully

[1] Blum H. A transformation for extracting new descriptors of shape. In: Dunn W, editor. Models for the perception of speech and visual form. Cambridge: MIT Press; 1967. p. 362–81. [2] Gelston S, Dutta D. Boundary surface recovery from skeleton curves and surfaces. Computer Aided Geometric Design 1995;12:27–51. [3] Yushkevich P, et al. Continuous medial representations for geometric object modeling in 2D and 3D. Image and Vision Computing 2003;21: 17–27. [4] Che WJ, Yang XN, Wang GZ. A dynamic approach to skeletonization. Journal of Software 2003;14(4):818–23 [in Chinese]. [5] Degen WLF. Exploiting curvatures to compute the medial axis for domains with smooth boundary. Computer Aided Geometric Design 2004;21:641–60. [6] Ramamurthy R, Farouki RT. 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. [7] Ramamurthy R, Farouki RT. 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. [8] Farouki RT, 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. [9] Choi HI, Han ChY, Moon HP, Wee N-S. New algorithm for medial axis transform of plane domain. Graphical Models and Image Processing 1997;59(6):463–83. [10] Teixeira RC. Medial axes and mean curvature motion I: Regular points. Journal of Visual Communication and Image Representation 2002;13: 135–55. [11] August AT, Zucker SW. On the evolution of the skeleton. Technical report DCS/TR 1179. Department of Computer Science, Yale University. [12] Ramanathan M, Gurumoorthy B. Constructing medial axis transform of planar domains with curved boundaries. Computer-Aided Design 2002; 35:619–32. [13] Sasaki S. Differential geometry. Tokyo: Kyoritsu Press; 1956 [in Japanese]. [14] Wang XC, Yu Y. An efficient surface–surface intersection algorithm based on geometry characteristics. Journal of Material Processing Technology 2002;123:191–6. [15] Li XY, Jiang H, Chen S, Wang XC. An efficient surface–surface intersection algorithm based on geometry characteristics. Computers & Graphics 2004;28:527–37.