- Email: [email protected]
Planar osculating arc splines
COMPUTER AIDED GEOMETRIC DESIGN ELSEVIER
Computer Aided Geometric Design 13 (1996) 653-671
Planar osculating arc splines ~ D.S. Meek*, D.J. Walton Department qf Computer Science, University (~fManitoba, Winnipeg, MB R3T 2N2, Canada Received August 1994; revised October 1995
Abstract A new method for drawing planar osculating arc splines, G 1 c u r v e s made by joining circular arcs, is described. The arc splines interpolate, match unit tangents, and match curvatures at the interpolation points. Arc splines produced by this new method match more general data, and thus give a more controlled curve, than arc splines produced from biarcs. The accuracy of arc splines produced by the new method is analyzed and compared to the accuracy of arc splines produced from biarcs when both are approximating smooth curves. When both have the same number of circular arcs, the new method produces an arc spline with a smaller error than the arc spline produced from biarcs. Keywords: Osculating arc spline; Accuracy of approximation
I. Introduction This paper studies the use of triarcs to produce planar osculating arc splines. An arc spline is a unit tangent continuous curve made by joining arcs. The arc spline is a desirable path for a numerically controlled cutting machine. Its advantages have been mentioned in recent papers (see for example (Meek, 1993)). Osculating arc splines will be derived from triarcs, and their accuracy of approximation will be compared to the accuracy of approximation of interpolating arc splines produced from biarcs. A biarc is a pair of circular arcs joined with a continuous unit tangent. A biarc can join two points and match the unit tangents at those two points. If a set of points and unit tangents is given, the biarcs joining adjacent pairs will produce an arc spline. The * The authors acknowledgethe financialsupport of the Natural Sciences and EngineeringResearch Council of Canada for this research. * Corresponding author. E-mail: [email protected]. 0167-8396/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0 167-8396(95)00053-4
654
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
biarc has one free parameter, and several ways have been proposed to set that parameter (see (Su, 1989; Meek, 1992)). A triarc is a set of three circular arcs joined with a continuous unit tangent. A triarc can join a pair of points matching unit tangents and matching curvatures at the two points. Not all curves are desirable for machining and the angle subtended by each arc is restricted in this paper. At interpolation points, the adjacent circular arcs from neighbouring triarcs pass through the same points, have the same tangent, and have the same curvature, so they are part of the same circle. Consequently, an arc spline produced by the triarc method with a large number of interpolation points has approximately two arcs per interpolation point. An arc spline produced by the biarc method also has approximately two arcs per interpolation point. As will be shown in Section 2, the triarc has one free parameter. Limits on the parameter are deduced, but no natural method for selecting the parameter in all cases has been discovered yet. An asymptotically optimal method for selecting the parameter in the important spiral case will be deduced in Section 3. A smooth curve will be taken to mean a curve with continuous third derivatives with respect to arc length. If an arc spline that approximates a given smooth curve is desired, the required unit tangents and curvatures can be calculated from the curve. The asymptotic error of an arc spline approximation produced by joining biarcs was analyzed in (Meek, 1995a). The corresponding asymptotic error for the arc spline approximation produced by joining triarcs is analyzed in Section 3. The asymptotic error in the arc spline produced from triarcs is smaller than that produced from biarcs when the two methods are used to approximate the same curve segment with approximately the same number of circular arcs.
2. Finding a triarc Given are points A, /3, unit tangents tA, tB, and the nonzero signed finite radii of the circles of curvature at A and B, RA, RB (see Fig. 1). Below only "general" types of triarcs are considered in detail. Many special cases are mentioned as they arise, but
JB
JA
CA
GB
Fig. 1. A triarc.
D.S. Meek, D.J. Walwn/ Computer Aided Geometric Design 13 (1996) 653~571
655
are not treated here. For example, the use of finite radii above prevents the first and third arcs from being straight lines. The convention on the sign of the radius is that the arc length is always positive, so the radius and the angle subtended by any given arc must have the same sign (counterclockwise rotation is positive). The centres of first and third circles, C A and C B , are CA = A + R A n A ,
(2.1)
Gu = B + Runu,
where n A and n g are the unit normals at A and B . Let JA be the point where the first and second arcs join and let J B be the point where the second and third arcs join. Let M be the centre of the circle from which the second arc is taken. To match tangents at JA and J g , M must be at the intersection of the lines J A V A and J B C B . Omit the special case where these two lines are parallel; this means the middle arc cannot be a straight line or a semicircle. M J A and M J B must be the same length, which is the radius of the second circular arc. Ignore the special case where CA and C B are identical and switch the labels A and B if necessary so that [RAI ~> IRBI. Let C = [ICB - CAll and assume a coordinate system is chosen so that CA :
~
,
CB = ~
.
The diagram can be reflected across the X - a x i s if necessary to make RA > 0. Notice that the given data determines circles through A and B . The subsequent analysis is based on the relation of these circles to each other. Some terminology follows. Let r be the signed radius of the middle arc. Let OA and 0B be the angles subtended by the first and third arcs. The angle 0 A will be used as the free parameter in each case. The following domains will be used to restrict the angles subtended by the arcs in the triarc. Define D1 to be the triangular region of points (x, y) such that 0 < x < 7r, 0 < y < 7r, x t- y < zr; D2 to be ( x , y ) such that -~r < :c < O, 0 < y < 7r, y - x < 7r; D3 to be ( x , y ) such that -Tr < x < O, -Tr < y < O, x + y > --~: and D4 to be (x, y) such that 0 < x < 7r, -Tr < y < 0, y - x > -Tr (see Fig. 2).
(O,x)
(-~,0)
U
(x,O)
(0,-~) Fig. 2. Domains DI, D2, D3, and D4.
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
656
(O,n)
(-rt,0)~
(n,0)
(O,-n)
Fig. 3. Plot of d(x, y) = p for various p. Let domain D be the square in the x - y plane bounded by the four lines y = - x + 7r, y = x + 7r, y = - x - 7r, and y = x - 7r. For ( x , y ) in D define siny-sinx_ d(x,y)-
y+x
cos 2 cos
(2.2)
A plot of d(x, y) = p for various p appears in Fig. 3. Lemma 2.1. For (x, y) in D and x nonzero, the curve d(x, y) = p, p > O, is equivalent to the curve
y = 2tan -I \l l- -p-(~ p cot 2 ) '
(2.3)
Proof. The equation d(x, y) : p can be written y+x COS
y-x -- p COS -
-
2 2 since cos(y - x ) / 2 ~ 0 in D. Expanding the cosines and collecting terms, tan y - 1 - P c o t 2 , 2 l+p from which Eq. (2.3) follows.
[]
Corollary. For (x, g) in D and g nonzero, the curve d(x, y) = p, p > O, is equivalent to the curve 1-p
x = 2 tan-I ( 1 - ~ p cot ~ ) .
(2.4)
D.S, Meek, D.J. Walton/ Computer Aided Geometric Design 13 (1996) 653~571
657
(O,n)
o/- \
\l/
/
\=.~
(-n,O)
(x,o)
\/o . # \ \ q / /I ./%'~' ,,,v, (O,-n) Fig. 4. Plot of e(x, y) = p for various p. For (x, y) in D define sin v+~ sin u - ~ ' 2 A plot of e(x, y) = p for various p appears in Fig. 4.
e(z, y) -
sin y + sin z sin(y - z)
(2.5)
L e m m a 2.2. lf p < 1 < q, point (a,b), a ~: b, a ¢ O, b ¢ O, in D such that d(a, b) = p and e(a, b) = q has components a=i2tan
-t
~/(1 - p ) ( q - 1) ~-~p)(q~l)'
and
b = i 2 t a n -1
~/(1 - p)(q + 1) (-1Tp)(q--1)"
Proof. Using the method of Lemma 2.1 on (2.5) for x ¢ y, the curve e(x, y) = q is equivalent to the curve
y=2tan-'(qq
+1 tan 2 ) . -1
Equating (2.3) and the above expression gives the result for a. The result for b can be obtained similarly. [] A detailed look at some general types of triarcs follows. The two cases of /~A and RB being the same sign, and of RA and RB being opposite signs will be examined.
Case 1: RA and R ~ are the same sign. By the previous discussion, both RA and RB are positive, and RA >1 RB. Ignore the special case RA = RB. Further, ignoring circles
D.£ Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
658
that are tangent to each other, there are three possibilities: the circle of curvature at A encloses the circle of curvature at B ; the two circles of curvature are separate and not enclosing; and the two circles of curvature intersect at two points. Let c~ and/3 be the angles in (-Tr, 7r) from the positive X-axis to the vectors A - CA and B - CB. In Case 1, both radii RA and RB are positive, so OA and OB are positive. Let x = o~ + Oa be the angle of JA -- C A and let y --/3 - 0B be the angle of J B -- CB. The angles x and y will be restricted somewhat arbitrarily so that (x, y) belongs to Dr, O2, D3, or D4. This avoids circular arcs with large angles and some special cases. The angle subtended by the second arc is y - x (see Fig. 1). Assume x ~ y as x = y would correspond to the special case of the second arc being a straight line. Finally, let R A -- R B
P -
C
RA + RB
'
q -
C
(2.6)
Note that both p and q are positive and p < q. The centre of the second arc M is the intersection of the two parametric lines
L l ( t ) = ~-
+
At M , t and u are sin y tC sin(y - x) but a t M ,
t=RA-r,
sin y - sin x sin(y - x)
ksinx/
and
u-
and
Lz(u) = ~
+
ksiny/"
sin x C, sin(y - x)
andu=Rs-r,
sot-u=RA-RB,
(2.8) and
RA -- RB C
'
or
d ( x , y) = d(o: q- OA, /3 -- OB ) = p,
(2.9)
where d(x, y) is defined by (2.2). OA is the free parameter of the triarc and p can be calculated from the given data, so (2.9) is an equation for OB in terms of OA. To find the radius of the middle arc, note that t + u = RA + RB -- 2r = Ce(x, y), where e(x, y) is defined by (2.5). It follows that r --
RA + RB 2
C
C
2 e(x,y)---- - ~ ( q - e ( x , y ) ) .
(2.10)
Equation (2.10) expresses the radius of the middle arc r in terms of OA. Consequently, Eqs. (2.9) and (2.10) give the whole triarc in terms of OA. There are three subcases labelled (a), (b), and (c) that arise from the relation between the two circles of curvature at A and B in Case 1.
Case l(a). The circle of curvature at A completely encloses the circle of curvature at B , or C < RA -- RB; from (2.6), I < p < q. The possible values of (x, y) are along the curve d(x, y) = p and, since p > 1, (x, y) is i n / ) 2 or 194 (see Fig. 3). In 0 2 and 04, e(x,y) is in the interval ( - 1 , 1) (see Fig, 4). Formula (2.10) shows that r is positive in D2 and 0 4 but the angle 9 - x subtended by the second arc is negative in/94. Thus, (x, y) cannot be in/94.
D.S. Meek, D.J. Walton/ Computer Aided Geometric Design 13 (1996) 653~571
659
(O,n)
(-n,o)
(0,0)
Fig. 5. Case l(a): Possible values for (x, y) in D2. If d(c~,/3) > p, then there are (x, y) such that d(x, y) = p and 0 A and 0B are positive (shown by the heavier part of the curve in Fig. 5). The lower bound on OA is 0. As OA increases, OB decreases, so the upper bound on OA occurs at 0B = 0. Using (2.4), the upper bound on on OA OA is upper bound is
OM = - c ~ + 2 t a n - '
-1 - p p cot ~ )
(2.1 l)
D
It is not known what value of OA would be best to use, but one possible choice would be the average of the limits, Oa = OM/2. In the error analysis of the asymptotic case, Section 3, a different and "optimal" OA is found. The range for r can be found from (2.101): C r < ? - ( q + l) :
RA + R s + C RA + RB + RA -- RB < = RA, 2 2
C r > ~(q-
RA + RB - C RA + RB 2 >
and 1) =
-- ( R A
-- R B )
:
I~B,
2
so RA > r > RB. The triarc in this case is a spiral. The circles of curvature at A and B suggest a spiral as there is a theorem due to Vogt (Guggenheimer, 1963) that circles of curvature of a spiral enclose all subsequent circles of curvature as one moves in the direction of increasing curvature. See Fig. 6 for two examples of spiral triarcs.
Case l(b). The circles of curvature at A and at B are separate and one does not enclose the other, RA + RB < C; from (2.6), p < q < 1. The possible values for (x, y) are along the curve d(x, y) = p and, since p < 1, (:r, y) is in D1 or D3 (see Fig. 3). In D1 with y > x, e ( x , y ) is in (1, oc) (see Fig. 4) and r is negative by (2.10); in D1 with y < x, e ( x , y ) is in ( - e c , - 1 ) and r is positive. In both cases, the angle y - x subtended by the second arc and its radius have opposite signs, so (x, y) is not allowed in Dl.
660
D.S. Meek, D.J. Walton/ ComputerAided Geometric Design 13 (1996) 653--671
CA
A Fig, 6. Case l(a): Spiral triarcs.
(-u,o)
(o,o) 0A
t
(o,-u) Fig. 7. Case l(b): Possible values for (x, y) in D3.
There are (x, y) such that d(x, y) = p and OA and 0B are positive (shown by the heavier part of the curve in Fig. 7) for any (a,/3) in D3. If d(a,/3) < p, then the lower bound on Oa is OM in (2.11) (see Fig. 7). If d(a,/3) >~ p, then the lower bound on OA is 0 (see Fig. 10). The upper bound on OA in both cases is - a . If y > x in D3, the angle subtended by the second arc is positive, and e(x, y) is in ( - o c , - 1 ) . Since q < 1, r is positive and the triarc forms a C-shaped curve. Further, from (2.10),
C RA + RB + C r> ~(q + l)= 2 > RA + RB.
D.S. Meek, D.J. Walton/Computer Aided Geometric Design 13 (1996) 653~571
661
Fig. 8. Case l(b): C-shaped and double S-shaped triarcs. (O,rO
(a,~) OB .--.i . ( x , y )
y=x
oA / (0,0)
(~,0)
Fig. 9. Case l(c): Possible values lor (x, y) in DI. If y < x in D3, the angle subtended by the second arc is negative, e(x, y) is in (1, oc), and the radius r is negative. The triarc is a double S-shaped curve. Further, from (2.10), C
r < ~-(q-
R A + R B -- C
1)=
2
See Fig. 8 for a C-shaped and two double S-shaped triarcs.
Case l(c). The circle of curvature at A intersects the circle of curvature at B at two places, R A -- R B < C < RA + RB; from (2.6) p < 1 < q. The possible values for (x, y) are along the curve d(x, y) = p and, since p < 1, (x, y) is in D1 or D3 (see Fig. 3). If y < x in D1, r is positive (as in Case l(b)) and the angle y - x subtended by the second arc is negative; this is not allowed. If y > x in D1 and e(x, y) > q, r is negative and y - x is positive; this is not allowed. The possible values of (x, y) in D1 are those along the curve d(x, y) = p for which y > x, and e(x, y) <~ q. Let (a, b) be the point in D1 such that d(a, b) = p and e(a, b) = q. L e m m a 2.2 gives the formula for finding this point. There are (x, y) such that d(x, y) = p and Oa and OB are positive (shown by the heavier part of the curve in Fig. 9) for any (a,/3) such that a < a and fl > b. If d(c~,fl) > p, then the lower bound on OA is OM in (2.11) (see Fig. 9). If d(c~,/3) ~< p, then the lower bound on OA is 0. The upper bound on OA is a--oL.
662
D.S. Meek, D.J. Walton/ ComputerAided GeometricDesign 13 (1996) 653~571 (-n,0)
(0,0)
(0,-re)
Fig. 10. Case l(c): Possible values for (x,y) in D3. CA
Ca
Fig. 11. Case l(c): C-shaped and double S-shaped triarcs. Since e(x, y) <~ q, r is positive and the triarc forms a C-shaped curve. Further,
C r=-~(q-e(x,y))
C RA ÷ RB -- C < - } - ( q - l) = 2
so r < RB < RA. This is a C-curve as shown in Fig. 1. The possible values for (x, y) in D3 are those along the curve d(x, y) = p with e ( x , y ) >1 q (see Fig. 10). Let (a,b) be the point in D3 such that d(a,b) = p and e(a, b) = q. L e m m a 2.2 gives the formula for finding this point. There are (x, y) such that d(x, y) = p and OA and 0B are positive (shown by the heavier part of the curve in Fig. 10) for any (c~,/3) such that c~ < a and/3 > b. If d(c~,/3) < p, then the lower bound on 0A is OM in (2.11). If d(c~,/3) ~> p, then the lower bound on OA is 0 (see Fig. 10). The upper bound on OA in both cases is a - o~. If y > x in D3, the angle subtended by the second arc is positive, and e(x, y) is in ( - o c , - 1 ) . Since q < 1, r is positive and the triarc forms a C-shaped curve. If y < x in D3, the angle subtended by the second arc is negative, e(x, y) is in (1, ~ ) , and the radius r is negative. The triarc is a double S-shaped curve. See Fig. I 1 for a C-shaped and two double S-shaped triarcs.
Case 2: RA and RB are opposite signs. By the above discussion, RA is positive and RB is negative; RA ~> --RB. The triarcs produced in this case are S-shaped.
D.S. Meek, D.J. Walton/ Computer Aided Geometric Design 13 (1996) 653-671
663
M
r-RA
JB
B RB
JA
Fig. 12. Case 2: S-shaped triarc. Let a be defined as before and/3 be the angle in (-Tr, 7r) from the positive X-axis to the vector CB - B (see Fig. 12). Let x = ct + OA be the angle of JA - CA as before and let y = / 3 + (--Os) be the angle of CB - JB. The angle subtended by the second arc is again y - x. The angles x and y will be restricted so that (x, y) belongs to D1, D2, D3, or D4, and assume x ¢ y as before. Recalling that RB is negative, let R A q- ( - - R B )
P-
C
R A -- ( - - R B )
'
q-
C
Note that p is positive, q is positive or 0, and q < p. The centre of the second arc M can be calculated as before using (2.8). Here t - u = RA + ( - R B ) , so the condition on (x, y) is again given by (2.9), and the radius of the second arc r is again given by (2.10). S-curves are possible when the circles of curvature at A and B are such that one encloses the other, the circles intersect, or the circles are separate and not enclosing. Only the case of circles being separate and not enclosing will be considered here; that is, assume q < p < 1. The possible values of (x, y) are along the curve d(x, y) = p and since p < 1, (x, y) is in Di or D3 (see Fig. 3). In D1 with y > x, e ( x , y ) is in (1,cx~) and r is negative; in Da with y < x, e ( x , y ) is in ( - ~ , - 1 ) and r is positive. In both cases, the angle y - x subtended by the second arc and its radius are opposite signs, so (x, y) is not allowed in DI. If (x, y) is in D3 and d(c~,/3) < p, there are (x, y) such that d(x, y) = p and OA and --OB are positive (shown by the heavier part of the curve in Fig. 13). The lower bound on 04 is 0, and OM in (2.11) is the upper bound. Since RA and RB are opposite in sign, the triarc produced is an S-shaped curve. See Fig. 14 for two S-shaped triarcs. There are many more cases to consider with triarcs than with biarcs. However, once the proper case has been identified, the amount of calculation for triarcs is similar to the amount of the calculation for biarcs. The calculation for triarcs proceeds as follows: R A and RB are given, so p and q can be calculated. Once the free parameter OA has been
664
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
(-rt,o)
(o,o)
d(x,y)= p ~ ~
,,4
/I (0,-Tt)
Fig. 13. Case 2: Possible values for (x, y) in D3.
B
CA
Fig. 14. Case 2: S-shaped triarcs.
chosen, OB can be calculated from (2.9) and Lemma 2.1. The radius r of the second circular arc can be calculated from (2.10), and its centre can be calculated from (2.7) and (2.8).
3. Asymptotic error for spiral triarcs A smooth curve generally has a finite number of curvature maxima and minima. If a smooth curve is subdivided finely, almost every segment will be a spiral, or a curve with monotone curvature. Thus, the asymptotic accuracy of the triarc approximation to a smooth curve depends on the accuracy of approximation of a triarc to a spiral. Some notation regarding spiral segments will now be established. Assume the spiral of positive increasing curvature from A to B , Q ( s ) , SA <<. s <<. sB, is parametrized in terms of arc length; the total arc length is h = S s -- 8a. The Frenet formulae for curves in the plane are t ' ( s ) = k ( s ) n ( s ) and n ' ( s ) = - k ( s ) t ( s ) . Using t = t ( s m ) for the unit
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
665
CA
R
R
JB
n
J
A=Q(sA)
t
Fig. 15. Triarc approximation to the spiral Q(s).
tangent vector at A = Q,(SA), n = n ( S A ) for the unit normal at A, the derivatives of Q ( s ) with respect to arc length are: Q ' ( s a ) = t,
Qtt(SA)
=
kn,
Qttt(SA)
=
-k2t + k'n,
where k = k(SA) > 0 and k' = ~ k ( S A ) > 0 are the curvature and its derivative with respect to arc length at SA. Later, k" is used to represent dd~sk(SA). It is convenient to use a coordinate system based on t and n as the X and Y axes and A as the origin (see Fig. 15). The radii of the circles of curvature at A and B a r e R A = 1 / k and R B = 1 / k ( s B ) , and the centres of the circles of curvature are CA = A + R A n and CB --= B + R B n ( s B ) (from (2.1)). The proofs of the following lemmas and theorem appear in more detail in (Meek, 1995b). Lemma
3.1. For the curve segment Q(s),
p = l + ~ - ~ k2 h2 + o(h
.
Proof. The vector CB - CA R A -- R B
B - A + R B n ( S B ) -- R A n RA -- RB
can be expanded by Taylor series to give C B -- C A
-kh +
( 2
12k k2
+ o(h 3) (3.J)
666
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653471
The norm of this vector is by (A.1)
c
: l-k2h2+o(h
),
(3.2)
24 and taking reciprocals of both sides of (3.2) gives RA -- RB
p-
R A C- R B
-
1 q-- k2 2 O(h3). -zTh + L t4
[]
Lemma 3.2. The angle a from C B - C A to A - CA, and the angle/3 from C B -- C A
to B - CB are
c~ = - 2
1--~-7h2 + O(h3)
and
/3 = kh2 +
12k'kk",]~h2 + O(h3).
Proof. The notation V x W shall be taken to mean the scalar V x W = IIvll. [[Wll sine = V x W u - VuWx, where ¢ is the angle from vector V = (W:) to vector W = (w:). The angles a and/3 can be found using (A.2) on sin(x-
C B -- C A
C
A -- C A __ × [[A--CAH
k h _ kk" hZ 2 ~ +O(h3)'
and s i n / 3 - C B - CA × = C lib - CBII ~h +
12k' h2 + O(h3)
where the unit vector (CB - C A ) / C is found by dividing (3.1) by (3.2), and the unit vector (13 - C B ) / I I B - CBlf is found by expanding - n ( s B ) . [] Lemma 3,2 shows that the total angle subtended by the three arcs,/3-c~, is kh + O(h2). Let OA = 01 kh be the angle in the first arc, let 02kh+O(h 2) be the angle in the second arc, and let 03kh + O(h 2) be the angle in the third arc. The three angles are not independent as 02 and 03 can be expressed in terms of 01, but subsequent expressions are made more compact by retaining 02 and 03. The following lemma shows how to calculate 03 from 01; 0z can be found from the relation 0 1 + 0 2 + 0 3 = 1. Lemma 3.3. Let OA :
Olkh, then the relation d(a + O A , f l - OB) = p implies that OB = 03kh + O(h 2) where 03--
1 301 3(I - 201)'
Proof. Using Lemma 3.1 for the value for p, Lemma 2.1, (A.4), and (A.5), k / 3 - - 0 B - - 6(1--201) h + O(h2)" "
(3.3)
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 6534571
667
Thus, 0B is 1 - 301 kh + O(h 0t3 = . 3(1 - 2 0 ~ )
= O kh +
[]
(3.4)
Since OA must be positive, 01 > 0 for small h. OB must also be positive, so (3.4) shows that 0 < 01 < 1/3. As 01 varies from 0 to 1/3, 03 varies from 1/3 to 0, and 02 varies from 2/3 down to a minimum of , f 3 / 3 at 01 = (3 - , / 5 ) / 6 and back up to 2/3 at 01 = 1/3. L e m m a 3.4. The radial difference between a spiral and an osculating circle is --Tgk' :~+ O(g4), where 9 is the arc length along the spiral from the point of contact. Proof. Let A be the point of contact. Using Taylor series expansions, k2
3 q_ O / 4~
g - -~-g-
O(sA q-g) = A +
k
~,9 )
(3.5)
'
gg2 + ~g3 + o(g4)
The centre of the circle of curvature Ca = A + ~n, so
IIQ(SA q - g ) - C A l l - -
1
k
k'
_~93 + 0 ( 9 4 ) .
The radial difference between the spiral and circular arc is k'
IIQ(sA + g) -- CAll -- RA = ----~93 + O(94).
[]
Arc length is related to angle by k ds =-- dO, so at the joins of the arcs in a tnarc, 9 = 01 h + O(h 2) and 9 = (1 - 03)h + O(h2). The radial difference between an arc and the curve it approximates will be called the radial error. 3.1. The radial error in the triarc approximation of a spiral of increasing curvature Q(SA + 9) is
Theorem
•
kt
--~-g3 + O(g4) (1 - 203)U
-4-~2
ifO~g<.O,h+O(h2), k'
(g _ Olh)2h _ _~g3 + O((g, h) 4)
E(9) = if
gU( h -
g)3
(3.6)
< g <. (1 - O3> +
_~ O((g, h) 4) if
(1 - 03)h + O(h 2) <~g <~ h + O(h2),
where the first, second, and third lines give the errors in the first, second, and third arcs. Proof. Error in the first arc. In the first arc, 9 varies from 0 to 01 h + O(h2). L e m m a 3.4 gives the asymptotic error in the first arc. The magnitude of the error is monotone
668
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
increasing with 9. The largest error in magnitude in the first arc occurs at g = 01 h+O(h 2) 1 btA3h3 and is -g,~ vl,o + O(h4). Error in the third arc. For the third arc, g varies from (1 - 0 3 ) h + O(h 2) to h + O(h2). Lemma 3.4 gives an asymptotic expression that shows that the error is monotone and its largest magnitude occurs at 9 = (1 - 03)h + O(h2). The largest error is lk'O]h3 + O(h4). Error in the second arc. In the second arc, 01h + O(h 2) ~< g ~ (1 - 03)h + O(h2). The radial error is
]IQ(sA + g ) - M I I - ~ =
IlQ(sa + g) - M{{ 2 - r 2 (3.7)
IIQ(SA + g) - M[I + r '
where M is (see Fig. 15) M=A+~
(--sinOA) = A +
1 (sin0A) 1--COS0A
+r
COS0A
1 ((1--rk)sinOA) k
"
1-- (1-- rk) cos 0A
The formula for r given in (2.10) can be expanded to
1 r -- k
(I-2o3)k' h 20zk 2
+ O(h2).
(3.8)
Asymptotic expansions give
M=A+
I/k+O(h)'
Q(sA+g)=a+
1
O(92)
,
1
IlO(sa + 9) - MII = ~ + O((9, h)l), and r = ~ + O(h), so the denominator of the expression on the right in (3.7) is 2
[[Q(sm + g ) - MII + r = ~ + O ( ( g , h ) ' ) , where the notation O((g, h) '~) stands for a power series in g and h starting with nth order terms. The radial error in (3.7) becomes
HQ,(SA + 9 ) - M l l - r =
(k+o((g,h)l))(llQ(sA +9)-Ml[2-r2).
(3.9)
The expression IIQ(SA + g) - M[]2 - r: can be written
;
g3+O(g4)+2(l_rk)
( ( Q ( s a + g ) _ a ) . ( - - s i n O a ~ + 1-- coSOA ) \ cos0a ,/
,~
"
The substitution OA = O1kh in the above and use of (3.5) yields k' (1 -- 203)k' ~.~g3 + "20--~ (g-Olh)2h+O((g'h)4)' where the expansion for 1 - rk from (3.8) was used. Now from (3.9) the error is
[IQ(*A + 9) - MII - ," ----- (1 - 402 203)k' (9 - O,h):h - gg3 k' ..}_O((g, h)4).
[]
D.S. Meek, D.Z Walwn/ ComputerAided GeometricDesign 13 (1996) 653--671 First Arc
Second Arc
669
Third Arc
Fig. 16. Radial error E(9 ) when approximating a short spiral. Some comments on the error formulae of Theorem 3.1 follows. The radial error in the first arc is negative, so the first circular arc is to the right of the spiral. The radial error in the third arc is positive, so the third circular arc is to the left of the spiral. Consequently, the second arc must cross the spiral. The maximum error of the first arc is at its end point JA and the maximum error of the third arc is at its end point ,/B. Since the error E ( 9 ) is continuous, the maximum error of the whole triarc can be found by looking at the error of the second arc. The error of the second arc has an internal maximum and an internal minimum. An example of a typical error is shown in Fig. 16. The labelling of the axes has been omitted intentionally because the example from which it was produced is only an approximation to the asymptotic case.
Theorem 3.2. The maximum of the leading term of the error in arc splines generated by the triarc method is minimized for OA = I ( 3 -- v/3)kh and the minimized maximum error is
1 (2q-3)3/2kth3+O(h4).
(3.10)
24 Proof. From the above discussion, the maximum error occurs in the second arc. Theorem A1 in (Meek, 1995b) shows that the maximum error in the second arc is minimized for small h when 01 = 03, as might be expected from symmetry. Eq. (3.3) with 01 = 03, gives 01 = ~ (3 - x/3). Thus, the maximum error is minimized when OA = ~ (3 -- x/~)kh. Substitution of this optimal 01 value into the error in the second arc in (3.6) gives the result in the theorem statement. [] In magnitude, this error from the triarc method is about 82% of the error from the biarc method. The biarc approximation has a maximum error of 3-~Uh3+O(h 4) (Meek, 1995a), which is approximately .00309 Uh 3 + O(h4). The maximum error in the triarc approximation (3.10) is approximately .00254 Uh 3 + O(h4).
4. Example of the asymptotic error An example spiral curve is a quadratic B6zier curve with control points (0,200), (150, 0), (300,0). This spiral was approximated by an arc spline formed from triarcs using the optimal OA from Theorem 3.2 and by an arc spline formed from Sabin's biarcs
670
D.XMeeL D.Z Walton/Computer Aided GeometricDesign 13 (1996) 653--671
Table 1 Accuracy of the arc spline approximation to the quadratic B6zier curve number of subdivisions
error in triarc arc spline
error in biarc arc spline
8
7.661 x 10-3
9.035
X 10 - 3
16
8.367 x 10-4
1.052
X 10 - 3
32
9.811 x 10-5
1.235
× 10 - 4
64
1.190 x 10-5
1.482 x 10-5
128
1.468
1.810
256
1.823 x 10-7
x 10 - 6
X 10 - 6
2.235 x 10-7
(Sabin, 1977) as described in (Meek, 1995a). Table I gives the number of times the curve was subdivided at the point of the mid-parameter value, The distance from the arc spline to the quadratic B6zier is estimated by taking the maximum of one hundred radial distances for both types of arc splines. The triarc produced arc splines are more accurate than the biarc produced arc splines. The error in the triarc produced arc spline is about 82% of the error in the biarc produced arc spline, and the accuracy of the triarc appears to be O(h 3) as predicted by Theorem 3.2. Note that although the arc spline produced from triarcs matches given curvatures, it is not a curvature continuous curve.
5. Conclusion Triarcs give a way to form arc splines that interpolate given points, match unit tangents at those points, and match curvatures at those points. This gives the designer good control over the arc spline. Triarcs have one free parameter. An optimal choice for that parameter was found for approximation of a short spiral. In cases where a smooth curve such as a B6zier curve is being approximated by an arc spline, it is easy to calculate the extra information (curvatures), that the triarc method requires. The calculation effort to find the arcs of a triarc is roughly the same as the effort to find the arcs of a biarc. Finally, in the limit the triarc gives a more accurate approximation to the smooth curve than the biarc for the same number of circular arcs.
Acknowledgment The authors appreciate the comments of two anonymous referees. They helped improve the presentation considerably.
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653~571
671
Appendix The following Taylor series expansions were used: +z=
l+~z+O(z2),
(A.I)
sin - I x = x + O(x 3),
(A.2)
sinx = z + O(x3),
(A.3)
t a n z : z + O(z3),
(A.4)
tan -1 z = z + O(z3).
(A.5)
References Guggenheimer, H.W. (1963), Differential Geometry, McGraw-Hill, Inc., New York. Meek, D.S. and Walton, D.J. (1992), Approximation of discrete data by G 1 arc splines, ComputerAided Design 24, 301-306. Meek, D.S. and Walton, DJ. (1993), Approximating quadratic NURBS curves by arc splines, Computer-Aided Design 25, 371-376. Meek, D.S. and Walton, D.J. (1995a), Approximating smooth planar curves by arc splines, J. Computational and Applied Mathematics 59, 221-232. Meek, D.S. and Walton, D.J. (1995b), Planar osculating arc splines, Technical Report 95/10, Department of Computer Science, University of Manitoba. Sabin, M.A. (1977), The use of piecewise forms for the numerical representation of shape, Technical Report 77/60, Computer and Automation Institute, Hungarian Academy of Sciences. Su, B.-Q. and Liu, D.-Y. (1989), Computational Geometry--Curve and Surface Modeling, G.-Z. Chang, Trans., Academic Press Inc., San Diego.
Computer Aided Geometric Design 13 (1996) 653-671
Planar osculating arc splines ~ D.S. Meek*, D.J. Walton Department qf Computer Science, University (~fManitoba, Winnipeg, MB R3T 2N2, Canada Received August 1994; revised October 1995
Abstract A new method for drawing planar osculating arc splines, G 1 c u r v e s made by joining circular arcs, is described. The arc splines interpolate, match unit tangents, and match curvatures at the interpolation points. Arc splines produced by this new method match more general data, and thus give a more controlled curve, than arc splines produced from biarcs. The accuracy of arc splines produced by the new method is analyzed and compared to the accuracy of arc splines produced from biarcs when both are approximating smooth curves. When both have the same number of circular arcs, the new method produces an arc spline with a smaller error than the arc spline produced from biarcs. Keywords: Osculating arc spline; Accuracy of approximation
I. Introduction This paper studies the use of triarcs to produce planar osculating arc splines. An arc spline is a unit tangent continuous curve made by joining arcs. The arc spline is a desirable path for a numerically controlled cutting machine. Its advantages have been mentioned in recent papers (see for example (Meek, 1993)). Osculating arc splines will be derived from triarcs, and their accuracy of approximation will be compared to the accuracy of approximation of interpolating arc splines produced from biarcs. A biarc is a pair of circular arcs joined with a continuous unit tangent. A biarc can join two points and match the unit tangents at those two points. If a set of points and unit tangents is given, the biarcs joining adjacent pairs will produce an arc spline. The * The authors acknowledgethe financialsupport of the Natural Sciences and EngineeringResearch Council of Canada for this research. * Corresponding author. E-mail: [email protected]. 0167-8396/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0 167-8396(95)00053-4
654
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
biarc has one free parameter, and several ways have been proposed to set that parameter (see (Su, 1989; Meek, 1992)). A triarc is a set of three circular arcs joined with a continuous unit tangent. A triarc can join a pair of points matching unit tangents and matching curvatures at the two points. Not all curves are desirable for machining and the angle subtended by each arc is restricted in this paper. At interpolation points, the adjacent circular arcs from neighbouring triarcs pass through the same points, have the same tangent, and have the same curvature, so they are part of the same circle. Consequently, an arc spline produced by the triarc method with a large number of interpolation points has approximately two arcs per interpolation point. An arc spline produced by the biarc method also has approximately two arcs per interpolation point. As will be shown in Section 2, the triarc has one free parameter. Limits on the parameter are deduced, but no natural method for selecting the parameter in all cases has been discovered yet. An asymptotically optimal method for selecting the parameter in the important spiral case will be deduced in Section 3. A smooth curve will be taken to mean a curve with continuous third derivatives with respect to arc length. If an arc spline that approximates a given smooth curve is desired, the required unit tangents and curvatures can be calculated from the curve. The asymptotic error of an arc spline approximation produced by joining biarcs was analyzed in (Meek, 1995a). The corresponding asymptotic error for the arc spline approximation produced by joining triarcs is analyzed in Section 3. The asymptotic error in the arc spline produced from triarcs is smaller than that produced from biarcs when the two methods are used to approximate the same curve segment with approximately the same number of circular arcs.
2. Finding a triarc Given are points A, /3, unit tangents tA, tB, and the nonzero signed finite radii of the circles of curvature at A and B, RA, RB (see Fig. 1). Below only "general" types of triarcs are considered in detail. Many special cases are mentioned as they arise, but
JB
JA
CA
GB
Fig. 1. A triarc.
D.S. Meek, D.J. Walwn/ Computer Aided Geometric Design 13 (1996) 653~571
655
are not treated here. For example, the use of finite radii above prevents the first and third arcs from being straight lines. The convention on the sign of the radius is that the arc length is always positive, so the radius and the angle subtended by any given arc must have the same sign (counterclockwise rotation is positive). The centres of first and third circles, C A and C B , are CA = A + R A n A ,
(2.1)
Gu = B + Runu,
where n A and n g are the unit normals at A and B . Let JA be the point where the first and second arcs join and let J B be the point where the second and third arcs join. Let M be the centre of the circle from which the second arc is taken. To match tangents at JA and J g , M must be at the intersection of the lines J A V A and J B C B . Omit the special case where these two lines are parallel; this means the middle arc cannot be a straight line or a semicircle. M J A and M J B must be the same length, which is the radius of the second circular arc. Ignore the special case where CA and C B are identical and switch the labels A and B if necessary so that [RAI ~> IRBI. Let C = [ICB - CAll and assume a coordinate system is chosen so that CA :
~
,
CB = ~
.
The diagram can be reflected across the X - a x i s if necessary to make RA > 0. Notice that the given data determines circles through A and B . The subsequent analysis is based on the relation of these circles to each other. Some terminology follows. Let r be the signed radius of the middle arc. Let OA and 0B be the angles subtended by the first and third arcs. The angle 0 A will be used as the free parameter in each case. The following domains will be used to restrict the angles subtended by the arcs in the triarc. Define D1 to be the triangular region of points (x, y) such that 0 < x < 7r, 0 < y < 7r, x t- y < zr; D2 to be ( x , y ) such that -~r < :c < O, 0 < y < 7r, y - x < 7r; D3 to be ( x , y ) such that -Tr < x < O, -Tr < y < O, x + y > --~: and D4 to be (x, y) such that 0 < x < 7r, -Tr < y < 0, y - x > -Tr (see Fig. 2).
(O,x)
(-~,0)
U
(x,O)
(0,-~) Fig. 2. Domains DI, D2, D3, and D4.
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
656
(O,n)
(-rt,0)~
(n,0)
(O,-n)
Fig. 3. Plot of d(x, y) = p for various p. Let domain D be the square in the x - y plane bounded by the four lines y = - x + 7r, y = x + 7r, y = - x - 7r, and y = x - 7r. For ( x , y ) in D define siny-sinx_ d(x,y)-
y+x
cos 2 cos
(2.2)
A plot of d(x, y) = p for various p appears in Fig. 3. Lemma 2.1. For (x, y) in D and x nonzero, the curve d(x, y) = p, p > O, is equivalent to the curve
y = 2tan -I \l l- -p-(~ p cot 2 ) '
(2.3)
Proof. The equation d(x, y) : p can be written y+x COS
y-x -- p COS -
-
2 2 since cos(y - x ) / 2 ~ 0 in D. Expanding the cosines and collecting terms, tan y - 1 - P c o t 2 , 2 l+p from which Eq. (2.3) follows.
[]
Corollary. For (x, g) in D and g nonzero, the curve d(x, y) = p, p > O, is equivalent to the curve 1-p
x = 2 tan-I ( 1 - ~ p cot ~ ) .
(2.4)
D.S, Meek, D.J. Walton/ Computer Aided Geometric Design 13 (1996) 653~571
657
(O,n)
o/- \
\l/
/
\=.~
(-n,O)
(x,o)
\/o . # \ \ q / /I ./%'~' ,,,v, (O,-n) Fig. 4. Plot of e(x, y) = p for various p. For (x, y) in D define sin v+~ sin u - ~ ' 2 A plot of e(x, y) = p for various p appears in Fig. 4.
e(z, y) -
sin y + sin z sin(y - z)
(2.5)
L e m m a 2.2. lf p < 1 < q, point (a,b), a ~: b, a ¢ O, b ¢ O, in D such that d(a, b) = p and e(a, b) = q has components a=i2tan
-t
~/(1 - p ) ( q - 1) ~-~p)(q~l)'
and
b = i 2 t a n -1
~/(1 - p)(q + 1) (-1Tp)(q--1)"
Proof. Using the method of Lemma 2.1 on (2.5) for x ¢ y, the curve e(x, y) = q is equivalent to the curve
y=2tan-'(qq
+1 tan 2 ) . -1
Equating (2.3) and the above expression gives the result for a. The result for b can be obtained similarly. [] A detailed look at some general types of triarcs follows. The two cases of /~A and RB being the same sign, and of RA and RB being opposite signs will be examined.
Case 1: RA and R ~ are the same sign. By the previous discussion, both RA and RB are positive, and RA >1 RB. Ignore the special case RA = RB. Further, ignoring circles
D.£ Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
658
that are tangent to each other, there are three possibilities: the circle of curvature at A encloses the circle of curvature at B ; the two circles of curvature are separate and not enclosing; and the two circles of curvature intersect at two points. Let c~ and/3 be the angles in (-Tr, 7r) from the positive X-axis to the vectors A - CA and B - CB. In Case 1, both radii RA and RB are positive, so OA and OB are positive. Let x = o~ + Oa be the angle of JA -- C A and let y --/3 - 0B be the angle of J B -- CB. The angles x and y will be restricted somewhat arbitrarily so that (x, y) belongs to Dr, O2, D3, or D4. This avoids circular arcs with large angles and some special cases. The angle subtended by the second arc is y - x (see Fig. 1). Assume x ~ y as x = y would correspond to the special case of the second arc being a straight line. Finally, let R A -- R B
P -
C
RA + RB
'
q -
C
(2.6)
Note that both p and q are positive and p < q. The centre of the second arc M is the intersection of the two parametric lines
L l ( t ) = ~-
+
At M , t and u are sin y tC sin(y - x) but a t M ,
t=RA-r,
sin y - sin x sin(y - x)
ksinx/
and
u-
and
Lz(u) = ~
+
ksiny/"
sin x C, sin(y - x)
andu=Rs-r,
sot-u=RA-RB,
(2.8) and
RA -- RB C
'
or
d ( x , y) = d(o: q- OA, /3 -- OB ) = p,
(2.9)
where d(x, y) is defined by (2.2). OA is the free parameter of the triarc and p can be calculated from the given data, so (2.9) is an equation for OB in terms of OA. To find the radius of the middle arc, note that t + u = RA + RB -- 2r = Ce(x, y), where e(x, y) is defined by (2.5). It follows that r --
RA + RB 2
C
C
2 e(x,y)---- - ~ ( q - e ( x , y ) ) .
(2.10)
Equation (2.10) expresses the radius of the middle arc r in terms of OA. Consequently, Eqs. (2.9) and (2.10) give the whole triarc in terms of OA. There are three subcases labelled (a), (b), and (c) that arise from the relation between the two circles of curvature at A and B in Case 1.
Case l(a). The circle of curvature at A completely encloses the circle of curvature at B , or C < RA -- RB; from (2.6), I < p < q. The possible values of (x, y) are along the curve d(x, y) = p and, since p > 1, (x, y) is i n / ) 2 or 194 (see Fig. 3). In 0 2 and 04, e(x,y) is in the interval ( - 1 , 1) (see Fig, 4). Formula (2.10) shows that r is positive in D2 and 0 4 but the angle 9 - x subtended by the second arc is negative in/94. Thus, (x, y) cannot be in/94.
D.S. Meek, D.J. Walton/ Computer Aided Geometric Design 13 (1996) 653~571
659
(O,n)
(-n,o)
(0,0)
Fig. 5. Case l(a): Possible values for (x, y) in D2. If d(c~,/3) > p, then there are (x, y) such that d(x, y) = p and 0 A and 0B are positive (shown by the heavier part of the curve in Fig. 5). The lower bound on OA is 0. As OA increases, OB decreases, so the upper bound on OA occurs at 0B = 0. Using (2.4), the upper bound on on OA OA is upper bound is
OM = - c ~ + 2 t a n - '
-1 - p p cot ~ )
(2.1 l)
D
It is not known what value of OA would be best to use, but one possible choice would be the average of the limits, Oa = OM/2. In the error analysis of the asymptotic case, Section 3, a different and "optimal" OA is found. The range for r can be found from (2.101): C r < ? - ( q + l) :
RA + R s + C RA + RB + RA -- RB < = RA, 2 2
C r > ~(q-
RA + RB - C RA + RB 2 >
and 1) =
-- ( R A
-- R B )
:
I~B,
2
so RA > r > RB. The triarc in this case is a spiral. The circles of curvature at A and B suggest a spiral as there is a theorem due to Vogt (Guggenheimer, 1963) that circles of curvature of a spiral enclose all subsequent circles of curvature as one moves in the direction of increasing curvature. See Fig. 6 for two examples of spiral triarcs.
Case l(b). The circles of curvature at A and at B are separate and one does not enclose the other, RA + RB < C; from (2.6), p < q < 1. The possible values for (x, y) are along the curve d(x, y) = p and, since p < 1, (:r, y) is in D1 or D3 (see Fig. 3). In D1 with y > x, e ( x , y ) is in (1, oc) (see Fig. 4) and r is negative by (2.10); in D1 with y < x, e ( x , y ) is in ( - e c , - 1 ) and r is positive. In both cases, the angle y - x subtended by the second arc and its radius have opposite signs, so (x, y) is not allowed in Dl.
660
D.S. Meek, D.J. Walton/ ComputerAided Geometric Design 13 (1996) 653--671
CA
A Fig, 6. Case l(a): Spiral triarcs.
(-u,o)
(o,o) 0A
t
(o,-u) Fig. 7. Case l(b): Possible values for (x, y) in D3.
There are (x, y) such that d(x, y) = p and OA and 0B are positive (shown by the heavier part of the curve in Fig. 7) for any (a,/3) in D3. If d(a,/3) < p, then the lower bound on Oa is OM in (2.11) (see Fig. 7). If d(a,/3) >~ p, then the lower bound on OA is 0 (see Fig. 10). The upper bound on OA in both cases is - a . If y > x in D3, the angle subtended by the second arc is positive, and e(x, y) is in ( - o c , - 1 ) . Since q < 1, r is positive and the triarc forms a C-shaped curve. Further, from (2.10),
C RA + RB + C r> ~(q + l)= 2 > RA + RB.
D.S. Meek, D.J. Walton/Computer Aided Geometric Design 13 (1996) 653~571
661
Fig. 8. Case l(b): C-shaped and double S-shaped triarcs. (O,rO
(a,~) OB .--.i . ( x , y )
y=x
oA / (0,0)
(~,0)
Fig. 9. Case l(c): Possible values lor (x, y) in DI. If y < x in D3, the angle subtended by the second arc is negative, e(x, y) is in (1, oc), and the radius r is negative. The triarc is a double S-shaped curve. Further, from (2.10), C
r < ~-(q-
R A + R B -- C
1)=
2
See Fig. 8 for a C-shaped and two double S-shaped triarcs.
Case l(c). The circle of curvature at A intersects the circle of curvature at B at two places, R A -- R B < C < RA + RB; from (2.6) p < 1 < q. The possible values for (x, y) are along the curve d(x, y) = p and, since p < 1, (x, y) is in D1 or D3 (see Fig. 3). If y < x in D1, r is positive (as in Case l(b)) and the angle y - x subtended by the second arc is negative; this is not allowed. If y > x in D1 and e(x, y) > q, r is negative and y - x is positive; this is not allowed. The possible values of (x, y) in D1 are those along the curve d(x, y) = p for which y > x, and e(x, y) <~ q. Let (a, b) be the point in D1 such that d(a, b) = p and e(a, b) = q. L e m m a 2.2 gives the formula for finding this point. There are (x, y) such that d(x, y) = p and Oa and OB are positive (shown by the heavier part of the curve in Fig. 9) for any (a,/3) such that a < a and fl > b. If d(c~,fl) > p, then the lower bound on OA is OM in (2.11) (see Fig. 9). If d(c~,/3) ~< p, then the lower bound on OA is 0. The upper bound on OA is a--oL.
662
D.S. Meek, D.J. Walton/ ComputerAided GeometricDesign 13 (1996) 653~571 (-n,0)
(0,0)
(0,-re)
Fig. 10. Case l(c): Possible values for (x,y) in D3. CA
Ca
Fig. 11. Case l(c): C-shaped and double S-shaped triarcs. Since e(x, y) <~ q, r is positive and the triarc forms a C-shaped curve. Further,
C r=-~(q-e(x,y))
C RA ÷ RB -- C < - } - ( q - l) = 2
so r < RB < RA. This is a C-curve as shown in Fig. 1. The possible values for (x, y) in D3 are those along the curve d(x, y) = p with e ( x , y ) >1 q (see Fig. 10). Let (a,b) be the point in D3 such that d(a,b) = p and e(a, b) = q. L e m m a 2.2 gives the formula for finding this point. There are (x, y) such that d(x, y) = p and OA and 0B are positive (shown by the heavier part of the curve in Fig. 10) for any (c~,/3) such that c~ < a and/3 > b. If d(c~,/3) < p, then the lower bound on 0A is OM in (2.11). If d(c~,/3) ~> p, then the lower bound on OA is 0 (see Fig. 10). The upper bound on OA in both cases is a - o~. If y > x in D3, the angle subtended by the second arc is positive, and e(x, y) is in ( - o c , - 1 ) . Since q < 1, r is positive and the triarc forms a C-shaped curve. If y < x in D3, the angle subtended by the second arc is negative, e(x, y) is in (1, ~ ) , and the radius r is negative. The triarc is a double S-shaped curve. See Fig. I 1 for a C-shaped and two double S-shaped triarcs.
Case 2: RA and RB are opposite signs. By the above discussion, RA is positive and RB is negative; RA ~> --RB. The triarcs produced in this case are S-shaped.
D.S. Meek, D.J. Walton/ Computer Aided Geometric Design 13 (1996) 653-671
663
M
r-RA
JB
B RB
JA
Fig. 12. Case 2: S-shaped triarc. Let a be defined as before and/3 be the angle in (-Tr, 7r) from the positive X-axis to the vector CB - B (see Fig. 12). Let x = ct + OA be the angle of JA - CA as before and let y = / 3 + (--Os) be the angle of CB - JB. The angle subtended by the second arc is again y - x. The angles x and y will be restricted so that (x, y) belongs to D1, D2, D3, or D4, and assume x ¢ y as before. Recalling that RB is negative, let R A q- ( - - R B )
P-
C
R A -- ( - - R B )
'
q-
C
Note that p is positive, q is positive or 0, and q < p. The centre of the second arc M can be calculated as before using (2.8). Here t - u = RA + ( - R B ) , so the condition on (x, y) is again given by (2.9), and the radius of the second arc r is again given by (2.10). S-curves are possible when the circles of curvature at A and B are such that one encloses the other, the circles intersect, or the circles are separate and not enclosing. Only the case of circles being separate and not enclosing will be considered here; that is, assume q < p < 1. The possible values of (x, y) are along the curve d(x, y) = p and since p < 1, (x, y) is in Di or D3 (see Fig. 3). In D1 with y > x, e ( x , y ) is in (1,cx~) and r is negative; in Da with y < x, e ( x , y ) is in ( - ~ , - 1 ) and r is positive. In both cases, the angle y - x subtended by the second arc and its radius are opposite signs, so (x, y) is not allowed in DI. If (x, y) is in D3 and d(c~,/3) < p, there are (x, y) such that d(x, y) = p and OA and --OB are positive (shown by the heavier part of the curve in Fig. 13). The lower bound on 04 is 0, and OM in (2.11) is the upper bound. Since RA and RB are opposite in sign, the triarc produced is an S-shaped curve. See Fig. 14 for two S-shaped triarcs. There are many more cases to consider with triarcs than with biarcs. However, once the proper case has been identified, the amount of calculation for triarcs is similar to the amount of the calculation for biarcs. The calculation for triarcs proceeds as follows: R A and RB are given, so p and q can be calculated. Once the free parameter OA has been
664
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
(-rt,o)
(o,o)
d(x,y)= p ~ ~
,,4
/I (0,-Tt)
Fig. 13. Case 2: Possible values for (x, y) in D3.
B
CA
Fig. 14. Case 2: S-shaped triarcs.
chosen, OB can be calculated from (2.9) and Lemma 2.1. The radius r of the second circular arc can be calculated from (2.10), and its centre can be calculated from (2.7) and (2.8).
3. Asymptotic error for spiral triarcs A smooth curve generally has a finite number of curvature maxima and minima. If a smooth curve is subdivided finely, almost every segment will be a spiral, or a curve with monotone curvature. Thus, the asymptotic accuracy of the triarc approximation to a smooth curve depends on the accuracy of approximation of a triarc to a spiral. Some notation regarding spiral segments will now be established. Assume the spiral of positive increasing curvature from A to B , Q ( s ) , SA <<. s <<. sB, is parametrized in terms of arc length; the total arc length is h = S s -- 8a. The Frenet formulae for curves in the plane are t ' ( s ) = k ( s ) n ( s ) and n ' ( s ) = - k ( s ) t ( s ) . Using t = t ( s m ) for the unit
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
665
CA
R
R
JB
n
J
A=Q(sA)
t
Fig. 15. Triarc approximation to the spiral Q(s).
tangent vector at A = Q,(SA), n = n ( S A ) for the unit normal at A, the derivatives of Q ( s ) with respect to arc length are: Q ' ( s a ) = t,
Qtt(SA)
=
kn,
Qttt(SA)
=
-k2t + k'n,
where k = k(SA) > 0 and k' = ~ k ( S A ) > 0 are the curvature and its derivative with respect to arc length at SA. Later, k" is used to represent dd~sk(SA). It is convenient to use a coordinate system based on t and n as the X and Y axes and A as the origin (see Fig. 15). The radii of the circles of curvature at A and B a r e R A = 1 / k and R B = 1 / k ( s B ) , and the centres of the circles of curvature are CA = A + R A n and CB --= B + R B n ( s B ) (from (2.1)). The proofs of the following lemmas and theorem appear in more detail in (Meek, 1995b). Lemma
3.1. For the curve segment Q(s),
p = l + ~ - ~ k2 h2 + o(h
.
Proof. The vector CB - CA R A -- R B
B - A + R B n ( S B ) -- R A n RA -- RB
can be expanded by Taylor series to give C B -- C A
-kh +
( 2
12k k2
+ o(h 3) (3.J)
666
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653471
The norm of this vector is by (A.1)
c
: l-k2h2+o(h
),
(3.2)
24 and taking reciprocals of both sides of (3.2) gives RA -- RB
p-
R A C- R B
-
1 q-- k2 2 O(h3). -zTh + L t4
[]
Lemma 3.2. The angle a from C B - C A to A - CA, and the angle/3 from C B -- C A
to B - CB are
c~ = - 2
1--~-7h2 + O(h3)
and
/3 = kh2 +
12k'kk",]~h2 + O(h3).
Proof. The notation V x W shall be taken to mean the scalar V x W = IIvll. [[Wll sine = V x W u - VuWx, where ¢ is the angle from vector V = (W:) to vector W = (w:). The angles a and/3 can be found using (A.2) on sin(x-
C B -- C A
C
A -- C A __ × [[A--CAH
k h _ kk" hZ 2 ~ +O(h3)'
and s i n / 3 - C B - CA × = C lib - CBII ~h +
12k' h2 + O(h3)
where the unit vector (CB - C A ) / C is found by dividing (3.1) by (3.2), and the unit vector (13 - C B ) / I I B - CBlf is found by expanding - n ( s B ) . [] Lemma 3,2 shows that the total angle subtended by the three arcs,/3-c~, is kh + O(h2). Let OA = 01 kh be the angle in the first arc, let 02kh+O(h 2) be the angle in the second arc, and let 03kh + O(h 2) be the angle in the third arc. The three angles are not independent as 02 and 03 can be expressed in terms of 01, but subsequent expressions are made more compact by retaining 02 and 03. The following lemma shows how to calculate 03 from 01; 0z can be found from the relation 0 1 + 0 2 + 0 3 = 1. Lemma 3.3. Let OA :
Olkh, then the relation d(a + O A , f l - OB) = p implies that OB = 03kh + O(h 2) where 03--
1 301 3(I - 201)'
Proof. Using Lemma 3.1 for the value for p, Lemma 2.1, (A.4), and (A.5), k / 3 - - 0 B - - 6(1--201) h + O(h2)" "
(3.3)
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 6534571
667
Thus, 0B is 1 - 301 kh + O(h 0t3 = . 3(1 - 2 0 ~ )
= O kh +
[]
(3.4)
Since OA must be positive, 01 > 0 for small h. OB must also be positive, so (3.4) shows that 0 < 01 < 1/3. As 01 varies from 0 to 1/3, 03 varies from 1/3 to 0, and 02 varies from 2/3 down to a minimum of , f 3 / 3 at 01 = (3 - , / 5 ) / 6 and back up to 2/3 at 01 = 1/3. L e m m a 3.4. The radial difference between a spiral and an osculating circle is --Tgk' :~+ O(g4), where 9 is the arc length along the spiral from the point of contact. Proof. Let A be the point of contact. Using Taylor series expansions, k2
3 q_ O / 4~
g - -~-g-
O(sA q-g) = A +
k
~,9 )
(3.5)
'
gg2 + ~g3 + o(g4)
The centre of the circle of curvature Ca = A + ~n, so
IIQ(SA q - g ) - C A l l - -
1
k
k'
_~93 + 0 ( 9 4 ) .
The radial difference between the spiral and circular arc is k'
IIQ(sA + g) -- CAll -- RA = ----~93 + O(94).
[]
Arc length is related to angle by k ds =-- dO, so at the joins of the arcs in a tnarc, 9 = 01 h + O(h 2) and 9 = (1 - 03)h + O(h2). The radial difference between an arc and the curve it approximates will be called the radial error. 3.1. The radial error in the triarc approximation of a spiral of increasing curvature Q(SA + 9) is
Theorem
•
kt
--~-g3 + O(g4) (1 - 203)U
-4-~2
ifO~g<.O,h+O(h2), k'
(g _ Olh)2h _ _~g3 + O((g, h) 4)
E(9) = if
gU( h -
g)3
(3.6)
< g <. (1 - O3> +
_~ O((g, h) 4) if
(1 - 03)h + O(h 2) <~g <~ h + O(h2),
where the first, second, and third lines give the errors in the first, second, and third arcs. Proof. Error in the first arc. In the first arc, 9 varies from 0 to 01 h + O(h2). L e m m a 3.4 gives the asymptotic error in the first arc. The magnitude of the error is monotone
668
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653-671
increasing with 9. The largest error in magnitude in the first arc occurs at g = 01 h+O(h 2) 1 btA3h3 and is -g,~ vl,o + O(h4). Error in the third arc. For the third arc, g varies from (1 - 0 3 ) h + O(h 2) to h + O(h2). Lemma 3.4 gives an asymptotic expression that shows that the error is monotone and its largest magnitude occurs at 9 = (1 - 03)h + O(h2). The largest error is lk'O]h3 + O(h4). Error in the second arc. In the second arc, 01h + O(h 2) ~< g ~ (1 - 03)h + O(h2). The radial error is
]IQ(sA + g ) - M I I - ~ =
IlQ(sa + g) - M{{ 2 - r 2 (3.7)
IIQ(SA + g) - M[I + r '
where M is (see Fig. 15) M=A+~
(--sinOA) = A +
1 (sin0A) 1--COS0A
+r
COS0A
1 ((1--rk)sinOA) k
"
1-- (1-- rk) cos 0A
The formula for r given in (2.10) can be expanded to
1 r -- k
(I-2o3)k' h 20zk 2
+ O(h2).
(3.8)
Asymptotic expansions give
M=A+
I/k+O(h)'
Q(sA+g)=a+
1
O(92)
,
1
IlO(sa + 9) - MII = ~ + O((9, h)l), and r = ~ + O(h), so the denominator of the expression on the right in (3.7) is 2
[[Q(sm + g ) - MII + r = ~ + O ( ( g , h ) ' ) , where the notation O((g, h) '~) stands for a power series in g and h starting with nth order terms. The radial error in (3.7) becomes
HQ,(SA + 9 ) - M l l - r =
(k+o((g,h)l))(llQ(sA +9)-Ml[2-r2).
(3.9)
The expression IIQ(SA + g) - M[]2 - r: can be written
;
g3+O(g4)+2(l_rk)
( ( Q ( s a + g ) _ a ) . ( - - s i n O a ~ + 1-- coSOA ) \ cos0a ,/
,~
"
The substitution OA = O1kh in the above and use of (3.5) yields k' (1 -- 203)k' ~.~g3 + "20--~ (g-Olh)2h+O((g'h)4)' where the expansion for 1 - rk from (3.8) was used. Now from (3.9) the error is
[IQ(*A + 9) - MII - ," ----- (1 - 402 203)k' (9 - O,h):h - gg3 k' ..}_O((g, h)4).
[]
D.S. Meek, D.Z Walwn/ ComputerAided GeometricDesign 13 (1996) 653--671 First Arc
Second Arc
669
Third Arc
Fig. 16. Radial error E(9 ) when approximating a short spiral. Some comments on the error formulae of Theorem 3.1 follows. The radial error in the first arc is negative, so the first circular arc is to the right of the spiral. The radial error in the third arc is positive, so the third circular arc is to the left of the spiral. Consequently, the second arc must cross the spiral. The maximum error of the first arc is at its end point JA and the maximum error of the third arc is at its end point ,/B. Since the error E ( 9 ) is continuous, the maximum error of the whole triarc can be found by looking at the error of the second arc. The error of the second arc has an internal maximum and an internal minimum. An example of a typical error is shown in Fig. 16. The labelling of the axes has been omitted intentionally because the example from which it was produced is only an approximation to the asymptotic case.
Theorem 3.2. The maximum of the leading term of the error in arc splines generated by the triarc method is minimized for OA = I ( 3 -- v/3)kh and the minimized maximum error is
1 (2q-3)3/2kth3+O(h4).
(3.10)
24 Proof. From the above discussion, the maximum error occurs in the second arc. Theorem A1 in (Meek, 1995b) shows that the maximum error in the second arc is minimized for small h when 01 = 03, as might be expected from symmetry. Eq. (3.3) with 01 = 03, gives 01 = ~ (3 - x/3). Thus, the maximum error is minimized when OA = ~ (3 -- x/~)kh. Substitution of this optimal 01 value into the error in the second arc in (3.6) gives the result in the theorem statement. [] In magnitude, this error from the triarc method is about 82% of the error from the biarc method. The biarc approximation has a maximum error of 3-~Uh3+O(h 4) (Meek, 1995a), which is approximately .00309 Uh 3 + O(h4). The maximum error in the triarc approximation (3.10) is approximately .00254 Uh 3 + O(h4).
4. Example of the asymptotic error An example spiral curve is a quadratic B6zier curve with control points (0,200), (150, 0), (300,0). This spiral was approximated by an arc spline formed from triarcs using the optimal OA from Theorem 3.2 and by an arc spline formed from Sabin's biarcs
670
D.XMeeL D.Z Walton/Computer Aided GeometricDesign 13 (1996) 653--671
Table 1 Accuracy of the arc spline approximation to the quadratic B6zier curve number of subdivisions
error in triarc arc spline
error in biarc arc spline
8
7.661 x 10-3
9.035
X 10 - 3
16
8.367 x 10-4
1.052
X 10 - 3
32
9.811 x 10-5
1.235
× 10 - 4
64
1.190 x 10-5
1.482 x 10-5
128
1.468
1.810
256
1.823 x 10-7
x 10 - 6
X 10 - 6
2.235 x 10-7
(Sabin, 1977) as described in (Meek, 1995a). Table I gives the number of times the curve was subdivided at the point of the mid-parameter value, The distance from the arc spline to the quadratic B6zier is estimated by taking the maximum of one hundred radial distances for both types of arc splines. The triarc produced arc splines are more accurate than the biarc produced arc splines. The error in the triarc produced arc spline is about 82% of the error in the biarc produced arc spline, and the accuracy of the triarc appears to be O(h 3) as predicted by Theorem 3.2. Note that although the arc spline produced from triarcs matches given curvatures, it is not a curvature continuous curve.
5. Conclusion Triarcs give a way to form arc splines that interpolate given points, match unit tangents at those points, and match curvatures at those points. This gives the designer good control over the arc spline. Triarcs have one free parameter. An optimal choice for that parameter was found for approximation of a short spiral. In cases where a smooth curve such as a B6zier curve is being approximated by an arc spline, it is easy to calculate the extra information (curvatures), that the triarc method requires. The calculation effort to find the arcs of a triarc is roughly the same as the effort to find the arcs of a biarc. Finally, in the limit the triarc gives a more accurate approximation to the smooth curve than the biarc for the same number of circular arcs.
Acknowledgment The authors appreciate the comments of two anonymous referees. They helped improve the presentation considerably.
D.S. Meek, D.J. Walton / Computer Aided Geometric Design 13 (1996) 653~571
671
Appendix The following Taylor series expansions were used: +z=
l+~z+O(z2),
(A.I)
sin - I x = x + O(x 3),
(A.2)
sinx = z + O(x3),
(A.3)
t a n z : z + O(z3),
(A.4)
tan -1 z = z + O(z3).
(A.5)
References Guggenheimer, H.W. (1963), Differential Geometry, McGraw-Hill, Inc., New York. Meek, D.S. and Walton, D.J. (1992), Approximation of discrete data by G 1 arc splines, ComputerAided Design 24, 301-306. Meek, D.S. and Walton, DJ. (1993), Approximating quadratic NURBS curves by arc splines, Computer-Aided Design 25, 371-376. Meek, D.S. and Walton, D.J. (1995a), Approximating smooth planar curves by arc splines, J. Computational and Applied Mathematics 59, 221-232. Meek, D.S. and Walton, D.J. (1995b), Planar osculating arc splines, Technical Report 95/10, Department of Computer Science, University of Manitoba. Sabin, M.A. (1977), The use of piecewise forms for the numerical representation of shape, Technical Report 77/60, Computer and Automation Institute, Hungarian Academy of Sciences. Su, B.-Q. and Liu, D.-Y. (1989), Computational Geometry--Curve and Surface Modeling, G.-Z. Chang, Trans., Academic Press Inc., San Diego.