- Email: [email protected]
Geometric control of G2-cubic A-splines
COMPUTER AIDED GEOMETRIC DESIGN ELSEVIER
Computer Aided Geometric Design 15 (1998) 261-287
Geometric control of G2-cubic A-splines M a r c o P a l u s z n y a,,,
R i c h a r d R. P a t t e r s o n b,l
a Centro de Computaci6n Grdfica y Geometrfa Aplicada, Facultad de Ciencias, Universidad Central de Venezuela, Apartado 47809, Los Chaguaramos, Caracas 1041-A, Venezuela b Department of Mathematical Sciences, Indiana Universi~ - Purdue Universi~ Indianapolis, 402 N. Blackford Street, Indianapolis, IN, 46202-3216, USA Received August 1995; revised July 1997
Abstract A-splines are piecewise algebraic curves, constructed in a control polygon that is a sequence of triangles meeting at the vertices. The arc in a given triangle is a cubic segment that joins these vertices and interpolates also the slope and curvature at each endpoint. This arc in each triangle is controlled by three other shape handles: an interior point that is interpolated, the slope at the interior point, and a tension parameter. The individual arcs can never be singular and are always convex. An inflection point of an A-spline must be placed at a junction point. © 1998 Elsevier Science B.V.
1. Introduction A-splines are piecewise algebraic curves. Sederberg proposed their use in geometric modeling (Sederberg, 1984); Bajaj coined the name. The A-splines discussed in this paper are curvature continuous. This property is equivalent to G2; see (Degen, 1988; Hohmeyer and Barsky, 1989; Garrity and Warren, 1991). We restrict here to cubic curves, since these have sufficient degrees of freedom to form G 2 splines, but few enough that we can give a geometrically intuitive interpretation for each of them. Two applications for splines are free-form geometric design and Hermite interpolation. We will discuss both of these applications for cubic A-splines. When using an A-spline for free-form design one begins with a control polygon as in Fig. 1; a sequence of * Corresponding author. E-mail: [email protected]. J E-mail: [email protected]. 0167-8396/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PH S01 6 7 - 8 3 9 6 ( 9 7 ) 0 0 0 3 1 - 9
262
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287 P
'p 0
Fig. 1. An A-spline.
triangles Ti such that two adjacent triangles have collinear edges. The spline interpolates joints Pzi, is tangent there to each segment Pzi-lPzi+l and has prescribed curvatures at each junction point. Each of these curvatures can be set numerically or relative to a default value. Each segment of the spline is convex, but zero curvature is allowed at any of its joints. In particular this means that inflection points may be placed where required and cannot occur otherwise. For example in Fig. 1, the curvature at P4 must be zero. Each segment can be controlled independently, keeping fixed the curvatures at its endpoints. The tangency and curvature interpolation conditions at the joints use six of the nine available degrees of freedom. The remaining three are managed as geometric handles for each individual segment. These are: an interior interpolation point B, the tangent line at B within limits imposed by the curvatures at the ends, and a tension parameter r. When cubic A-splines are used for high accuracy Hermite interpolation, the points P2i are chosen on a given curve, and then using the tangent lines at the Pzi the control polygon is formed around it. If inflection points are present they must be selected to be among the P2i. The curvatures at the P2i are computed and become additional input, as well as another point B in each triangle and its tangent. Using the results of (Degen, 1993, 1995) one can see that for each fixed value of the tension parameter r the A-spline we construct here approximates the original curve to order 8. We start with an overview of some of the previous work on G 2 splines, beginning with parametric splines. A construction of G 2 conic splines is possible, but they can not interpolate curvatures given at the joints and have only one (global) shape handle (Pratt, 1985; Farin, 1989). Pottman describes G 2 splines created from elements, each of which is composed of two conic arcs pieced together (Pottmann, 1991b). These curves have local control, but can never have zero curvature and have fewer shape handles than cubic splines. Cubic spline interpolation using parametric cubics as described by Farin produces curves that interpolate points but not tangents or curvatures (Farin, 1993). The choice of a knot sequence is required, whose effect on the curve is not well enough understood. A global system of linear equations must be solved.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
263
The direct G 2 cubic splines ('y-splines) described in (Farin, 1993) interpolate only the endpoints and tangents. An alternative construction for these curves is given in (Pottmann, 1991a). A local scheme for the use of polynomial cubic curves for high accuracy (order 6) Hermite interpolation was proposed by (de Boor, H011ig and Sabin, 1987). The method involves solving a system of nonlinear equations. Shirman and S6quin extend this construction with more stress on the application to free-form design (Shirman and S6quin, 1992). They describe G 2 splines created from elements, each of which is composed of two cubic arcs pieced together. Degen extends the construction to rational curves and achieves interpolation of order 8 (Degen, 1993, 1995). The most difficult part concerns the existence and reality of the curves. Sederberg (1984) proposed the use of algebraic splines and Warren (1991) expanded on his ideas. Li, Hoschek and Hartmann (1990) gave a construction that produces high order geometric continuity at each junction point but uses only one of the available degrees of freedom for shape control. Levin and Nadler (1995) extended these ideas to give methods for interpolating points that lie on the boundary of a convex region with an oval of an algebraic curve. Bajaj and Xu (1992, 1996) discuss tangent, curvature and higher order geometric continuity of A-splines. They use all of the available degrees of freedom (subject to one inequality), not for direct shape control, but for higher order approximation at the junction points and for least squares approximation of additional data points within each triangle. They also show how algebraic splines can be constructed in a sequence of parallelograms or rectangles (Bajaj and Xu, 1995). The authors studied the construction of cubic algebraic splines in (Paluszny and Patterson, 1993, 1994a), using all but one of the possible degrees of freedom for direct shape control. Both tangent continuous and curvature continuous splines were constructed, but the construction of curves that were not necessarily globally convex was hindered by not using all the available degrees of freedom. In this paper all of the available degrees of freedom are used for direct shape control. We find a larger family than Bajaj and Xu were able to use, and the restriction to globally convex curves is removed. In fact we characterize the maximal family of cubic curve segments that satisfy the interpolation, tangency, and curvature conditions at the given points. There have been two main objections raised against using implicit splines. First it was assumed to be impossible to control the shape of a curve in a geometrically intuitive way by adjusting the coefficients. This paper should lay to rest that objection; we show how to use all of the available degrees of freedom for geometrically meaningful shape control. Second it was asserted that rendering and robustness problems would never be solved. This objection is not dealt with in this paper, but some progress has been made; see the references cited in Section 9. Naturally, parametric curves are easier to render and methods for creating them extend to space curves and tensor product surfaces. By contrast, implicit curves provide easy control of convexity and displaying a family of solutions certainly solves the problem of existence and reality of solutions (among which will be any parametrizable curves that exist). The conclusion is to provide a hybrid system
264
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
that takes advantage of the strengths of both parametric and implicit curves (Marsilio et al., 1994). The organization of the paper is as follows. In the next section individual cubic segments are studied and the interpolation point B is introduced. Section 3 describes the parameter space for controlling the segment and introduces the slope parameter m. In Section 4 a cubic discriminant A0 is defined. The graph of A0 = 0 in the parameter space delineates the portion of the parameter space that corresponds to useful curves. This acceptable region is studied further in Section 5. Section 6 studies the tension parameter r. In Section 7 we explore the use of G 2 cubic A-splines for free form geometric modeling and some proposed defaults for the initial settings of the shape parameters. Section 8 discusses high order Hermite interpolation and Section 9 contains a brief discussion of graphing implicit curves. Notation: This paper uses three different types of coordinate systems. They are distinguished notationally as follows: Triples P[s, t, u] in square brackets refer to the homogeneous coordinates of the point P, triples in round brackets P ( s , t, u) to the barycentric coordinates of P, and pairs P ( x , y) in round brackets to Euclidean coordinates.
2. Behavior of a cubic segment
In this section we investigate the behavior of a cubic curve segment in an individual triangle To. The segment is defined by an implicit equation with ten coefficients; however Sederberg observed that the interpolation of endpoints and end tangencies in each triangle causes four of the ten coefficients to be zero (Sederberg, 1984). The equation for the ith segment can then be expressed in terms of the barycentric coordinate system defined by the ith triangle in the control polygon as F = asZu + bsu 2 - cst 2 - dt2u + e s t u - f t 3 = O.
The coefficients can be associated to trisection points of the triangle as in Fig. 2. The remaining six coefficients can be used for geometric control of the curve segment. One can think of the graph of F = 0 as the zero'th level curve of a functional triangular surface patch, but the multinomial coefficients have been absorbed into a, b , . . . , f to simplify notation. Some conditions on the signs of these coefficients are necessary. It is simple to check that a and b must have the same sign in order for the curve not to cross the segment POP2. We must make sure that the curve enters the triangle at/90. If we choose Euclidean coordinates so that the vertices are P0(0, 0), PI (mj, 0) and P2(x2, y2) (see Fig. 3) and compute the branch expansion of F = 0 in a neighborhood of P0, we obtain C Y2
---x Y=ax~
2
+.--.
From this it is clear that c and a must have the same sign for the curve to enter the triangle at P0. A similar result holds for coefficients b and d. Thus it is necessary to have a, b, c, d of the same sign to have any hope of finding an appropriate arc inside the triangle from P0 to P2, so we now impose this condition.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
265
P0/ s=l
1.1=1
Fig. 2. Association of coefficients to points of the triangle.
..... "'"'"'"'"'"'"'"'""~.0I"//~/xl, x
Po(O, o) Fig. 3. Branch expansion of F ----0 near 19o.
If in addition f is not zero and if the sign of f is the same as that of c and d the curve will not cross either PoP1 or P1P2. In this case the only boundary points of the triangle that the curve passes through are Po and P2. The branch expansion shows that Y=
c Y2 X 2
ax-~
is an osculating parabola to the curve at the origin. Its curvature at the origin is easily computed to be
~ = 2 cy2 a x2"
Hence in the notation of Fig. 4 we have the following result. Proposition 1. If F = 0 is an algebraic curve of degree 3 in which the coefficients of s 3 and 82t a r e O, its curvature at Po is given by
2 c hi a 9o~
4c_A a g o~
where A is the area of triangle PoP1P2. Similarly if the coefficients of u 3 and tu 2 are O, its curvature at P2 is given by dho =4 d A N,1 = 2 ~ if-'}ff-~.
266
M. Paluszny, R.R. Patterson / ComputerAided Geometric Design 15 (1998) 261-287
-,.
os S
Fig. 4. Euclidean measurements of the triangle. It turns out that having coefficients a, b, c, d, f all of the same sign guarantees that the curve will have an arc from P0 to P2 and no inflection point or singular point inside the triangle (although it can have a singular point somewhere outside the triangle). The following result is due to Bajaj and Xu (1992, 1996). Theorem 2. For a, b, e, d, f all o f the same sign, any straight line that passes through 191 and crosses the line segment PoP2 intersects the curve one and only one time (counting multiplicities) in the interior o f the triangle PoP~ P2. Furthermore the curve segment can have no inflection point in the interior o f the triangle. The following result generalizes part of the above theorem by removing the condition on the coefficient f .
Proposition 3. l f a cubic curve contains a nonsingular arc that enters the triangle at Po tangent to PoP1 and leaves at P2 tangent to PIP2, then there are no inflection points on this arc in the interior o f the triangle.
Proof. A real cubic curve has either one or three inflection points. In the case of three, they are collinear (Salmon, 1879). Since the cubic segment in the triangle has the same concavity at P0 and P2, if it has any inflection points in the triangle it must have two. The line through them either separates P0 and P2 or it does not. If it does separate/90 and P2 then after the curve crosses the line twice (see Fig. 5a) it is trapped on the P0 side and cannot get to P2 without another inflection. If the line does not separate P0 and P2 (see Fig. 5b), the center portion of the segment, which must be concave up, cannot be joined to P0 and P2 without crossing the line again. D We now begin the analysis of how the coefficients can be used for geometric control of the curve. We can use two of the remaining degrees of freedom to specify curvatures at the endpoints. Referring to Proposition 1, define #0 = c / a and #l = d/b. These can be used to control the curvatures at the joints. From the branch expansion, both #0 and/zl must be nonnegative. The values e = #0a and d ~1 b can be substituted into F. =
267
M. Paluszny, R.R. Patterson / ComputerAided Geometric Design 15 (1998) 261-287
(a)
(b)
P2
Fig. 5. Inflection points inside the triangle.
There are three remaining degrees of freedom to be specified. For the first of these, a point B(s0, to, u0) in the triangle is selected. Requiring the curve to interpolate B allows one to solve and determine a specific value e0 for the coefficient e: eo = - ( a s 2 u o + bsou~ - a # o s o t ~ -
b#lt2uo - f t 3)/sotouo.
With the additional substitution of e0 for e, the equation of F becomes asZu + bsu 2 - a#ost 2 - b#lt2u + eostu - f t 3 = O,
(1)
where #0, Pl, so, to and uo are now geometrically meaningful parameters. It remains now to provide geometric meaning for the remaining parameters a, f and b, and to do so in such a way as to obtain the largest possible family of useful curve arcs.
3. The [a, f, b]-parameter
plane
Many useful curves exist for which the sign of f does not agree with that of a and b. If f = 0 the curve passes through P1, and for f of opposite sign to that of a, b one gets a branch of the cubic coming into the triangle. This is not a problem so long as there is a segment from Po to P2 that is traced, say from Po to P2, and the two curve branches do not touch (the curve remains nonsingular inside the triangle). Bajaj and Xu do not allow these extra curves because they are using least squares methods to find the best fitting curve segment to data points inside the triangle--a significantly different problem. With #o, #1 and B chosen, the cubic depends homogeneously on the parameters a, f and b. One can think of the projective plane of [a, f, b] as the parameter plane for the family of cubics of the form (1) with curvatures determined by #0 at P0 and #l at P2 and interpolating the point B. The shaded region in Fig. 6 consists of those points A(a, f, b) for which ab > 0. By Theorem 2, at least every point in the triangle P a P f P b gives a cubic that has a convex segment connecting P0 and P2. Using B for the first shape handle, the other two arise from the following observation.
Proposition 4.
Cubic segments with the same #0, ~1 and B that have the same tangent line at B correspond to parameter points A that lie on a straight line in parameter space.
268
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 2 6 1 - 2 8 7
Fig. 6. Preliminary domain of A.
P,
, ~ ' ~ [ s 1, O, Ul] Fig. 7. Interpolating cubic arc. Proof. The equation of the tangent line to F at B is
Fs(B)s + Ft(B)t + F,~(B)u = 0. This line meets the PoP2 line at the point Q[sl, 0, Ul], where [Sl,0,?~l] = [ - Fu(B),O, Fs(B)]; see Fig. 7. The point Q lies outside the segment PoP2 because the cubic arc can have no inflection points inside the triangle. Such points Q are parametrized by m ~
81 q- U l 81 - - u 1
for-l~
--
1. Thus
81 -Jr Ul __ - -
81 -- '/3,1
Fu(B) - F,(B) f u ( B ) + F~(B)'
which simplifies to
asg [m(#oto2 + u~) - (#oto2 - u2)] + ft~ [rn(so + uo) - (so - uo)] + bug[m(p,tg + s~) + (#it 2 - sg)] = 0.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
269
Fig. 8. Preliminary strategy for controlling cubic segments. This is the equation of a straight line in the [a, f, b]-plane.
[]
The line in Proposition 4 is called an m-line. It is remarkable, but easy to check, that all of the m-lines pass through a common point
I[t30uo(souo- #,t~),souo(#o#,t 4 - s~u~),t30so(souo- #ot~)] • The point I exists except in the case in which/-to and #, happen to be the curvatures of a conic arc that solves the problem. Indeed, if the notation 7 = t~/souo is introduced, the coordinates of I are / [~o7(1 - #17), #o#,72 - 1, u~7(1 - #o7)] ;
(2)
and all three coordinates of I are zero in the case #0 = #1 = 1/7. In this case the cubic F is reducible:
F = (Tsu - t 2) (as + f T t + bu). The quadratic factor is a solution to the problem: it satisfies the interpolation conditions at P0 and P2, has the curvatures determined by #0 and #1, and interpolates B. The parameters a, f and b have no affect on the conic. Proposition 4 suggests a strategy for controlling the cubic segments: select a point B to be interpolated, select a tangent line at B by means of choosing an m-line in the pencil through I, and then for the third shape handle use the remaining parameter r to move the parameter point A along this line. The situation is much like polar coordinates; see Fig. 8. However, the limits on m and r can be extended beyond the PaPyPb triangle to a natural boundary that is given by a certain curve Ao = 0 in the parameter plane.
4. The cubic discriminant /t
There is a universal polynomial A of degree twelve whose variables are the ten coefficients of a general homogeneous cubic polynomial in three variables (Salmon, 1879).
270
M. Paluszny, R.R. Patterson / Computer Aided Geometric
Design 15 (1998) 261-287
It has the property that when d > 0 the cubic has one circuit (as a curve in projective space), when L < 0 the cubic has two circuits, and when n = 0 the cubic is singular or reducible. In the case under consideration, in which four of the coefficients are zero, the formula for n simplifies (Paluszny and Patterson, 1993). After the substitutions c = pea and d = ~1 b, A reduces further to A=a3b3(-e3f3+27abf4+e4f2po-36abef3po+8abe2f2& + 16a2b2f2& + e4f2p1 - 36abef3p1 - e5f,uopl +46abe2f2popul - 8abe3f&l
- 24a2b2f2&I
- 24a2b2f2p&
- 16a*b*ef&~
- abe4p&: + 64a*b*ef&T
+ 16a2b2f2& - 16a2b2efp&
+ 8abe2f2p: - 8abe3fp& - 8a2b2e2&T - 16a3b3&f
- 8a2b2e2p&: + 32a3b3&:
- 16a3b3&;)
= a3b3A0. The additional substitution e = ea creates the sixth degree homogeneous polynomial Ao(a, f, b), which is too long to reproduce here, but can easily be found with a computer algebra program. The graph of A0 = 0 in the [a, f, b]-plane will help to delineate a larger region in which the parameter points A(u, f, b) that determine useful cubic arcs may be selected. Proposition 5. When not both po and ~1 equal l/y, I never lies inside the triangle P,PfPb.
the point I is a singularity
of
do.
Proof. The coordinates of I satisfy A0 = 0, (no), = 0, (A,), = 0, and ( Ao)b = 0, so it is a singular point of multiplicity at least two. Using (2) one sees it is impossible for q all three coordinates of I to have the same sign, so I lies outside the triangle. The graph of A0 differs considerably depending on the relative sizes of po, pi and l/y. To simplify notation, introduce the variables va = ,LLOY - 1 and v1 = ply - 1 by making the substitutions PO = (uo + l)SO210& in Ao. In this notation I are
and
PI = (211+ ~)souo/$,
the formula for A0 becomes shorter again and the coordinates
1 [&or4, -(vov1 + uo + u,,&L;,
so&].
of
(3)
Because y > 0 and ~0, ~1 3 0, one has ~0, vi > - 1. This region in (~0, vl)-space is the meta-parameter space for studying Ao. The origin corresponds to the case in which I does not exist. In order to determine the type of the singularity I for A,, the Hessian of A0 must be examined as a quadratic form (Patterson, 1988a):
HeNAo)
=
(no),,(l)
(Ao)&)
(Ao)ab(I)
(AoLd
(Ao>.~J(~
(Ao)fb(r)
(Ao)ab(l)
(Ao),b(r)
(Ao)bb(l)
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
271
vI
(0,o)
Vo
(-1,-1)
Fig. 9. Graph of h and the meta-parameter space.
The following computation can be carried out using computer algebra. First each entry of Hess(Ao) is found to have the factor Sot 8,12 0 u o8-2 n , where 2 21 . h ( v o , Vl) = v 2 - 2 v o v , + v~ + 2 v 2 v , + 2 v o v 2 + VoV
Assuming that vo ¢ 0 and h ¢ O, the remaining matrix of Ao can then be reduced by a sequence of identical row and column operations to the matrix
(ooo) 0 0
h 0
0 0
.
Therefore the graph of h (see Fig. 9) separates the meta-parameter space into the region in which h > 0 and I is a crunode and the shaded region in which h < 0 and I is an acnode (isolated singularity). This result holds on the line Vl = 0 since no assumption was made about vl. It holds also on the line v0 = 0 since a corresponding result can be obtained by reducing the matrix using the (3, 3)-element of the Hessian as the first pivot instead of the (1, 1)-element. On the graph of h = 0, I is a third order singularity.
272
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
(b)
Fig. 10. Example graphs of Ao when vo, Vl > O and (a) vo = ~Ul
=
3, (b) vo = 3, vj = 4.2.
Here are some examples of the graph of Ao. In case v0, vl > 0, I is a crunode in the region ab > 0. In Fig. 10a, vo = Vl = 3. In this special case of vo = vl the sextic Ao factors into a quartic and a double linear factor whose graph is a double tangent. The point I is circled. In Fig. 10b, v0 = 3 and Vl = 4.2. In these examples the graph of F = 0 interpolates B ( 1 / 4 , 1/2, 1/4). In case v0, vl < 0, I also lies in the region ab > 0. In Fig. l l a , v0 = vl = - 0 . 5 . Again Ao reduces to a quartic and a double line. In Fig. 1 lb, vo = - 0 . 4 5 and vl = - 0 . 5 . In these examples (v0, vl) lies in the region in which I is an acnode and the graph of F = 0 interpolates the barycenter B ( 1 / 3 , 1/3, 1/3). Notice that the graph of A0 = 0 sometimes enters the PaPfPb-triangle. A curve F = 0 that corresponds to a point A on the graph of Ao = 0 in the triangle is singular, but its singular point does not occur
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
273
•
f
"--.-,~.~..,,,..~
(b)
®
Fig. 11. Example graphs of Ao when vo, vl < 0 and (a) vo = vj -- - 0 . 5 , vl = - - 0 . 5 .
(b) vo = - 0 . 4 5 ,
274
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
Fig. 12. Example graph of Ao when vj < 0 < vo.
inside the PoPl P2-triangle. These singular curves F = 0 must be the polynomial and rational solutions found by (de Boor, Hollig and Sabin, 1987; Degen, 1993, 1995). Fig. 12 shows one example in which v0 = 1 > 0, v~ = - 0 . 5 < 0, and B is the barycenter. Now I is a crunode that does not lie in the region ab > O.
5. The acceptable region Parameter points A that produce useful segments of the corresponding curves can be chosen to lie in the region ab > 0 below the PaPyPb triangle and above the graph of A0 = 0. (The previous sentence needs to be understood in the sense of projective geometry, in which the region below the triangle is connected to the region above the triangle through the line at infinity.) This region is called the acceptable region. By
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
275
1)l
4 1)0
Fig. 13. Areas in the meta-parameter space.
sampling parameter points in other regions one finds that the corresponding curves do not satisfy the requirements. The shape of the acceptable region depends on the point (v0, Vl ) in the meta-parameter space. This point is determined by the curvatures at the endpoints as reflected in #0 and #l and on the interpolated point B. Number the areas in meta-parameter space as in Fig. 13. For the case of area 1, in which v0, Vl > 0, the acceptable region is four sided with two curved edges and I at a vertex, as in Fig. 14a or 14b. The two curved edges in general do curve outward as shown, so all of the points in the region are visible from the point I. Most of them lie in the acceptable quadrilateral determined by I and the three vertices Pa, Pf and Pb. For the case of area 2a or 2b, in which v0, vl < 0, an arc of A0 = 0 crosses the lines a = 0 and b = 0, separating the point I from the triangle. In this case the acceptable region is three sided as in Fig. 15a or 15b. The curved edge in general does arch away from the triangle. Most of the acceptable points lie in the acceptable triangle CaPfSb (Fig. 15a) or SaPfCb (Fig. 15b). The remaining cases in which one of v0, vl is positive and the other negative are similar to the previous, except that the curve forming the lower edge of the three sided region passes through Pa or Pb, as in Fig. 16. Now imagine the m-line through I sweeping across the acceptable region. In many cases it is necessary to find an intersection of this line with the graph of A0 = 0. By Bezout's Theorem there can be as many as four intersections; finding the correct one would be possible but slow. For this reason it is proposed to restrict to the acceptable quadrilateral or the acceptable triangle. This compromise misses some acceptable curves but allows many more than the restriction to the original triangle PaPfPb. In order to find the limits on m in each case it is necessary to compute the value of m at each vertex of the acceptable region. This is easiest for area 1. In the v0, '/31 notation, the equation of an m-line is
a.o~o [m(.o~o * ~o + .o) -- ("o~o -- ~o + "o)] + St~ [m(~o + "o) + ~o -- ~o] + b~o~ [m(~o~, + ~o + ~o) + (~o~, + ~o - ~o)] = o.
(4)
276
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
Fig. 14a. In area 1, to cover the acceptable quadrilateral, m varies from the value m a ~-
80Vo -- UO + 80 sOVO + uo + 80
for the m-line through Pa to m b .~
'/z0v I ~ u 0
-- 80
UoVl -{- UO ÷ sO
for the m - l i n e t h r o u g h Pb. T h e n r v a r i e s so as to m o v e A ( a , f , b) f r o m I to o n e o f s e g m e n t s P a P I o r P I P b , d e p e n d i n g o n m (see Fig. 14a, b). T h e c u t - o f f v a l u e is
mf-
8o
-
the
u0
so + uo
which is the value for the m-line through Pf. This verifies the information in the first row of the table in Fig. 17.
M. Paluszn); R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
277
Fig. 14b. Acceptable region when vo, vl > 0 and (a) I lies below the triangle, (b) I lies above the triangle.
For areas 2, 3 and 4, since some of the vertices of the acceptable region are intersections of the graph of A0 = 0 with the lines a = 0 and b = 0, it is necessary first to find the coordinates of the points Sa, Sb, Ca, C6 and determine which ones are actually the vertices. When a = 0, A0 factors
-- f 74 Vl ( f Tto + buov l ) 3 ( f~/to + buo + b~tovl ) ( fa/tovo + buo v l + buovo v l ). The cubic factor corresponds to the point Tb[0, uovl,--')'to], where A0 = 0 has the line a = 0 as an inflectional tangent; see the graphs of Ao = 0. The next factor determines a point Sb[O, uo(1 + Vl),--Tto], and the last factor determines Cb [0, uovl (1 + vo),--Ttovo]. Similarly when b = 0, Ao factors
- f @ v o ( f T t o + asovo) 3 (aso + f~/to + asovo)(asovo + fTtovl + asovovl).
278
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
St)
$a Cb
!
I
Fig. 15a. The corresponding points are Ta[-Tto, sovo, 0], Sa[-Tto, so(1 + vo), 0]
and
Ca[-TtOVl, sovo(1 + vL), 0].
In order to distinguish the points SD, Cb in the graph of Ao one can use the cross ratio. Given four distinct points P/J1,0],/95[0, 1], U[u, 1], V[v, 1] on the projective line, their cross ratio is given by
cR(Ps, Pb, u, v) =
1
0
0
1
u
1
v
1
v
0 u
1 1
1 v
0 1
u
The points U and V are separated by Pf, Pb if and only if C R < O. If U and V are not separated by Pf, Pb then they occur in the order P f P b U V if v / u > 1 and in the order P f P b V U if v / u < 1.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
279
,I
Cb
Sa
Fig. 15b. Acceptable region when vo, vt < 0, I is an acnode, and (a) I lies below the triangle, (b) I lies above the triangle.
Applied to the problem at hand one has 1
0
CR(Ps P~,S~,C~)= uo(1 + v~) -Tto '
0 uo(1 + v l )
1 -~to
0
1
uovl(1 + vo) -Ttovo 1
0
u0vl(l + vo)
-~/tovo
(vo+ 1)v~ (v, +l)vo Thus the four points P f , Pb, Sb, C6 occur in this order if and only if 'vt < v0; in area 2a. Similarly the points Py, Pa, So, Ca Occur in this order if and only if v0 < vl; in area 2b. It follows that the vertices of the acceptable triangle are P f , Co and Sb in area 2a, and p f , S,~ and Cb in area 2b.
280
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
cb
Fig. 16. Acceptable region when v~ < 0 < vo.
Area
1
Values of v0, v~
Vertices
vo, vl > 0
P~,Ps,Pb, I
Domain of m m a to m b
2a
vl < v o < 0
C~,Ps, Sb
Tt~ca to ft2sb
2b
vo < vl < 0
Sa, Ps, Cb
r o s a to T/Zcb
3
vl < 0 < v o
P,~, P~ , Sb
m y to m s b
4
vo<0
S,~,Ps,Pb
rns~ t o r a f
Fig. 17. Limits of the acceptable region. Two more rows of the table in Fig, 17 are completed by finding the values of m for which the m-line passes through Sa, Sb, Ca, Cb. These are respectively rosa = - 1 , m s b ~- 1,
281
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
(b) Fig. 18. Region in which I lies in Case 3.
mca
=
(s0
-
uo)(vo
-
-
uovo
l
(SO + U0)(V0 -- Vl) + U0V0U1
and
mcb =
(so
-
u0)(
o -
-
sovo
l
(SO + U0)(V0 -- Vl) -- 80V0Vl
.
In the case of area 3, v0 > 0, vl < 0 as in Fig. 16, so the vertices of the acceptable triangle are Pa, Py and either Sb or Cb. Since both components of Cb have the same sign, Cb lies between P f and Pb. However the signs of the components of Sb are opposite, so 8b is the comer of the acceptable triangle. In Case 3 the only other thing needed to complete the table is to verify that the point I is located in the [a, f, b]-projective plane in such a way that the extreme vertices of the acceptable triangle as seen from I are always P I and Sb, never Pa. Note that the line t3of + (1 + Vl)Sou2b = 0 through Pa and Sb and the line b = 0 determine two regions in the projective plane, and that the triangle P ~ P I P b lies in one of them. It is enough to check that for v0 > 0 and vl < 0, I is contained in the same region as the triangle. Figs. 18a and 18b illustrate this for two possible locations of the point Sb. The computation for area 4 is analogous to that for area 3.
6. The tension parameter r Once #0, #1, B and m have been chosen, there is a family of curves F = 0 determined by the remaining parameter r. These curves have three points of contact at each of Po and P2 since they have the same tangent line and curvature, and two points of contact at B since they have the same tangent line. Because all cubics through eight fixed points pass through a ninth fixed point (Walker, 1950), it follows that all curves in this family share another point J. If J does not lie on the curve arc from Po to P2 inside the triangle, then for any two curves of the family, one lies entirely to one side of the other except at the endpoints and B. And in fact, most of the time r serves as a tension parameter, pulling the curve towards or away from the polygon determined by the sides u = 0, s = 0 and the tangent line at B.
Proposition 6. /f/z0, #1, B are fixed so that (vo, vl ) lies in area 1 or 2, and m is within the bounds delimited by the acceptable region, then the ninth point o f intersection o f any two cubics does not lie on the curve arc inside the triangle.
282
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
Proof. For (v0, V,) in area 1 or 2 let Q(ao, 0, bo) be the point where the m-line intersects the line P~Pb. Parametrize the m-line by (1 - r ) l + rQ. Choose two distinct values rl and r2 and let F1 and F2 be the corresponding cubics. The t-resultant of F1 and F2 is
]~(r 1
-
-
T2)3(VoVl -'~ V0 -Jr-U1)S3U3(S/Z0
-
-
2tSO)2(aoa2s + b0/32u),
where the factor k depends only on B, and a and/3 in the last linear factor LI are
a = aosov~ + bouov~ + bouovov~, /3 = aosov 2 + bouov~ + aosov2v,. The factors s 3 and u 3 correspond to P2 and P0, the factor (suo - uso) 2 corresponds to B, and the linear factor gives a line L, = 0 through Pl and J. Similarly the s-resultant of Fl and F2 has a linear factor L2 that represents a line through Po and J :
L2 = -aosouo(vovl + vo + vl ) a t + to(aosovo - bouovj )flu. Solving Ll and L2 simultaneously one finds the coordinates of J
d [bosouo(vovl + vo + vt)/32, to(bouov, - aosovo)a/3,-aosouo(vovl + vo + Vl)a2]. When m is such that Q lies between Pa and Pb, ao and b0 are positive, so the first and third coordinates of J have opposite signs and J cannot even lie in the region su > O, thus cannot lie inside the triangle PoP1P2. In particular when (v0, vl) is in area 1 and A is in the acceptable quadrilateral J cannot lie inside the triangle. However when m is varied so that Q no longer lies between P~ and Pb, J can enter the triangle. This situation is easier to examine using actual values of a0 and b0, which are determined from (4):
ao = uo(m(uovl + uo + so) + uov, + uo - so) bo = - s o ( m ( s o v o
and
+ uo + ~o) - ~ovo + ~o - so).
The coordinates of J become J[j,~, jr, ju], where
js = - s o ( m ( s o v o + uo + so) - sovo + uo - so)/3 2, j~ = t o ( m ( ~ 0 + ~0) + ~0 - so ) a / 3 ,
j~
=
- u o ( m ( u o v , + uo + so) + uov, + uo
-
s o ) a 2,
and where now a -- m ( ( s o + ~ o ) ( ~ o
-~l)
- so~o~l) - ((so - - o ) ( ~ o
- ~,) - ~o~o~,)
and /3 -- m ( ( s o + uo)(vo - v , ) + u o v o v , ) - ((so - uo)(vo - ~,) - u o ~ o v , ) .
Thus one sees that as m sweeps the acceptable region, js = 0 only when m = inca, which makes/3 ----0, or m = m~, which makes the other factor zero. For (vo, Vl) in area 2a, race is a bound of the domain of m, while in area 2b, mc~ is outside the domain of m. So for area 2, /3 can vanish only at one end of the domain of m. For m = m~
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
283
the coordinates of J are [0, t0vo(1 + v l ) , u o ( v o v i + vo + Vl)]. For (vo, vl) in area 2, jt and ju have the same sign, so J lies on the P1P2-side of the triangle. Thus J can enter the Pi P2-side of the triangle only on the other branch of the cubic. A similar argument holds for the PoPl-side of the triangle. So for (v0, vl) in area 2, J cannot lie on the arc from P0 to P2. The only case that remains to be proved is when (v0,vl) lies in area 1 and A is in the acceptable region but outside the acceptable quadrilateral. This proof is omitted because of the decision not to use these curves anyway. [] So r acts as a tension parameter whenever both p0,,ttl > 1/7 ((U0,'UI) lies in area 1) or #0, #1 < 1/7 ((v0, vl) lies in area 2). For some choices of (v0, vl) in areas 3 and 4, J does lie on the curve arc from P0 to P2 inside the triangle and can even be B. In practice in these cases it is observed that motion of A along such an m-line is essentially doing very little to change the graph of F = 0 anyway. Here are some illustrations of the use of the tension parameter r. For v0, vl > 0 (area 1), when A approaches one of the bottom curved boundaries of the acceptable region, the curve develops a comer as in Fig. 19a. In the extreme case of A actually on the bottom boundary the curve has a crunode in triangle PoP1P2; when A = I the crunode is at B. For vo, Vl < 0 (area 2), when A approaches the bottom curved boundary of the acceptable region the curve fits tightly against the tangent line at B as in Fig. 19b.
7. G 2 cubic A-splines in free-form m o d e l i n g
In this section we discuss some defaults that can be used so that as soon as a designer has constructed a control polygon as in Fig. 1, a curve will appear. First consider curvature. Curvature is an invariant of Euclidean geometry and introducing it is somewhat unfortunate since all the other constructions with F are being done in affine geometry. Without knowing the scale of a figure it is impossible to appreciate the significance of a number as curvature or radius of curvature at a point. One consequence of having defaults for the curvatures is that a slider can then be used to control the curvature at a junction point. In order to construct default curvatures at the junction points, consider two adjacent triangles T i - i and Ti and the graph of su - t 2 = 0 in each. In the notation of Fig. 4, the curvature of the segment in Ti-l at P2i is 2hzi_z/g2i_l while that of the curve in Ti is 2h2i+1/92i. Ordinarily these curvatures will be different. However if we define )ki--
h2i-2g~i 2 h2i+lg2i_ 1
and use V/~iisu - t 2 = 0
in Ti-1
and
- - 1s u
_ t2 = 0
in Ti,
284
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
(b)
Fig. 19. Example cubic arcs for A near the boundary of the acceptable region. these conics both have the same curvature at P2~, namely the geometric mean of the original two curvatures ~i = 2 ~ / h 2 i - 2 h 2 i + l g2i-lg2i
Although we cannot generally create a G 2 c o n i c spline in this way, we can use ~i as the default curvature at P2i for a cubic spline. If the spline is open, at the endpoints we can use the curvature of s u - t 2 = 0; hi
~O=2g-~
2
and
~'~=
h2n-2
g22n_l"
The next defaults are the interior points B~ that are interpolated; they can be chosen to be the barycenters. Choosing B~ also has the effect of selecting an acceptable region
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
285
in parameter space for control of the ith segment of the spline. The limits on the slope parameter mi at B~ are determined by the extreme vertices of the acceptable region as seen from I, whose position depends on Bi and the curvatures at the endpoints of the ith segment; see Fig. 17. The parameter mi can also be controlled by a slider, with the default being an intermediate value. The tension ri can then be controlled by yet another slider that moves the point A on the segment of the m-line in the projective plane. The limits on ri are determined by the value of mi and the line segments that bound the acceptable quadrilateral or triangle, and like m~ it can default to an intermediate value.
8. G2 cubic A-splines for Hermite approximation Suppose a smooth curve C is given, either parametrized or at least in such a way that it is possible to select points on it and compute tangents and curvatures. Choose points P2i on C, including any inflection points. The tangents at the P2i determine a sequence of triangles. Choose another point Bi on C in each triangle. We have seen how to construct cubic segments from P2i through B to P2i+2 that have the same tangent and curvature at each end, the same tangent at Bi and one remaining degree of freedom r. Let Mi be the midpoint of P2i and P2i+2. Because both C and the cubic segment (for any fixed r) are convex, rays from M~ provide a transverse direction along which distance can be measured between the two curves. Degen has shown that the maximum of this distance goes to zero as the eighth power of the distance h from PEi to P2i+2 as the partition is refined (Degen, 1993, 1995). (Actually he measures distance h in a parameter space, but this is equivalent.) This approximation order could theoretically be increased to nine if a parameter could be found to replace r and increase the order of contact at B or if the curve were to interpolate two interior points instead of just one. If a curve segment comes from a conic it will be reproduced exactly; cubic A-splines have the property of quadratic precision. Because of the endpoint and tangency conditions on the conic, in its equation the coefficients of s 2, '/.Z2, 8 t and tu will be zero. We can then normalize the equation to the form 7su - t 2 = 0. Because B lies on the curve, = t~/(souo). The curvatures of this conic at Po and P2 are given by n0 -
2 hi
"y go
and
/~1 :
2 ho
---
gF
So /~0 :
~1 :
l/"y.
This is the case in which I does not exist and the cubic is reducible; see Section 3. The conic is one of the factors and is the arc that is traced from P0.
9. Graphing implicit curves There are a number of methods available for graphing implicit curves: (Bajaj et al., 1988; Chandler, 1988; Patterson, 1988b; Sederberg et al., 1989; Dobkin et al.,
286
M. Paluszn3, R.R. Patterson/ Computer Aided Geometric Design 15 (1998) 261-287
1990; Bajaj and Xu, 1994; Marsilio et al., 1994; Paluszny and Patterson, 1994b; Taubin, 1994). The simpler of the two methods described by Taubin was used to graph h and A0 in Figs. 9, 10-12, 14-16. Cubic A-spline segments can be traced, since their behavior inside each triangle is so predictable. The fastest and most robust method that we have found to graph them uses a ray that emanates from the midpoint of the PoP2 segment and rotates in steps from P0 to P2. The ray is intersected with the tangent line to the curve at the previously computed point, and then a Newton correction is made to find the next point on the curve. This method was used for Fig. 19.
Acknowledgments The authors wish to thank W.L.E Degen for reading a draft of this paper and providing many helpful suggestions. This research was partially supported by the NSF/Conicit U.S.-Latin America Cooperative Science Program, grant INT-9304208.
References Bajaj, C., Hoffmann, C., Hopcroft, J. and Lynch, R. (1988), Tracing surface intersections, Computer Aided Geometric Design 5, 285-307. Bajaj, C. and Xu, G. (1992), A-splines: local interpolation and approximation using Ck-continuous piecewise real algebraic curves, Computer Science Technical Report, CAPO-92-44, Purdue University. Bajaj, C. and Xu, G. (1994), Rational spline approximations of real algebraic curves and surfaces, in: Dikshit, H.P. and Micchelli, C.A., eds., Advances in Computational Mathematics, C, World Scientific Publishing Co, Approximations and Decomposition Series, in press. Bajaj, C. and Xu, G. (1995), Regular algebraic curve segments (III) - Applications in data fitting, manuscript, Purdue University. Bajaj, C. and Xu, G. (1996), Data fitting with cubic A-splines, in: Gigante, M. and Kunii, T., eds., Proceedings of Computer Graphics International, CGI'94, Melbourne, Australia, World Scientific Publishing Co. Chandler, R. (1988), A tracking algorithm for implicitly defined curves, IEEE Computer Graphics and Applications 8, 83-89. de Boor, C., HOllig, K. and Sabin, M. (1987), High accuracy geometric Hermite interpolation, Computer Aided Geometric Design 4, 269-278. Degen, W. (1988), Some remarks on Brzier curves, Computer Aided Geometric Design 5, 259-268. Degen, W. (1993), High accurate rational approximation of parametric curves, Computer Aided Geometric Design 10, 293-313. Degen, W. (1995), High accuracy approximation of parametric curves, in: Dahlen, M., Lyche, T. and Schumaker, L.L., eds., Mathematical Methods in CAGD II1, to appear. Dobkin, D., Levy, S., Thurston, W. and Wilks, A. (1990), Contour tracing by piecewise linear approximations, ACM Transactions on Graphics 9, 389--423. Farin, G. (1989), Curvature continuity and offsets for piecewise conics, ACM Transactions on Graphics 8, 89-99. Farin, G. (1993) Curves and Surfaces for Computer Aided Geometric Design, Academic Press, New York.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
287
Foley, T., Goodman, T. and Unsworth, K. (1994), An algorithm for shape preserving parametric interpolating curves with G 2 continuity, in: Laurent, P.-J., Le Mrhautr, A. and Schumaker, L.L., eds., Curves and Surfaces in Geometric Design, A K Peters, Wellesley, MA, 249-259. Garrity, T. and Warren, J., (1991), Geometric continuity, Computer Aided Geometric Design 8, 51-65. Hohmeyer, M. and Barsky, B. (1989), Rational continuity: parametric, geometric and Frenet frame continuity of rational curves, ACM Transactions on Graphics 10, 335-359. Levin, D. and Nadler, E. (1995), Convexity preserving interpolation by algebraic curves and surfaces, Numerical Algorithms 9, 113-139. Li, J., Hoschek, J. and Hartmann, E. (1990), G n ~-functional splines for interpolation and approximation of curves, surfaces and solids, Computer Aided Geometric Design 7, 209-220. Marsilio, R, Paluszny, M., Patterson, R. and Rodriguez, O. (1994), Ambiente h~rido de disefio gr~ifico: splines algebraicos y paramrtricos, in: Memorias XX Conferencia Latinoamericana de Informdtica, Mexico City, Mexico, 147-156. Paluszny, M. and Patterson, R. (1993), A family of tangent continuous cubic algebraic splines, ACM Transactions on Graphics 12, 209-232. Paluszny, M. and Patterson, R. (1994a), A family of curvature continuous cubic algebraic splines, Technical Report #94-8, Department of Mathematical Sciences, Indiana University - Purdue University Indianapolis. Paluszny, M. and Patterson, R. (1994b), G2-continuous cubic algebraic splines and their efficient display, in: Laurent, P.-J., Le Mrhautr, A. and Schumaker, L.L., eds., Curves and Surfaces in Geometric Design, A K Peters, Wellesley, MA, 353-359. Patterson, R. (1988a), Parametric cubics as algebraic curves, Computer Aided Geometric Design 5, 139-159, Patterson, R. (1988b), Parametrizing and graphing nonsingular cubic curves, Computer-Aided Design 20, 615~i23. Pottmann, H. (1991a), A projective algorithm for curvature continuous rational splines, in: Farin, G., ed., NURBS for Curve and Surface Design, SIAM, New York. Pottmann, H. (1991b), Locally controllable conic splines with curvature continuity, ACM Transactions on Graphics 10, 366-377. Pratt, V. (1985), Techniques for conic splines, in: Proceedings of SIGGRAPH 85, ACM, New York, 151-159. Salmon, G. (1879), Higher Plane Curves, Hodges, Foster and Figgis, Dublin. Sederberg, T. (1984), Planar piecewise algebraic curves, Computer Aided Geometric Design 1, 241-255. Sederberg, T., Zhao, J. and Zundel, A. (1989), Approximate parametrization of algebraic curves, in: Strasser, W. and Seidel, H.-P., eds., Theory and Practice of Geometric Modeling, Springer, 33-54. Shirman, L.A. and S6quin, C.H. (1992), Procedural interpolation with curvature-continuous cubic splines, Computer-Aided Design 24, 278-286. Taubin, G. (1994), Rasterizing algebraic curves and surfaces, IEEE Computer Graphics and Applications 14, 14-23. Walker, R. (1950), Algebraic Curves, Dover, New York. Warren, J. (1991), Several notions of geometric continuity for implicit plane curves, Monograffas de la Academia de Ciencias de Zaragoza 2, 121-129.
Computer Aided Geometric Design 15 (1998) 261-287
Geometric control of G2-cubic A-splines M a r c o P a l u s z n y a,,,
R i c h a r d R. P a t t e r s o n b,l
a Centro de Computaci6n Grdfica y Geometrfa Aplicada, Facultad de Ciencias, Universidad Central de Venezuela, Apartado 47809, Los Chaguaramos, Caracas 1041-A, Venezuela b Department of Mathematical Sciences, Indiana Universi~ - Purdue Universi~ Indianapolis, 402 N. Blackford Street, Indianapolis, IN, 46202-3216, USA Received August 1995; revised July 1997
Abstract A-splines are piecewise algebraic curves, constructed in a control polygon that is a sequence of triangles meeting at the vertices. The arc in a given triangle is a cubic segment that joins these vertices and interpolates also the slope and curvature at each endpoint. This arc in each triangle is controlled by three other shape handles: an interior point that is interpolated, the slope at the interior point, and a tension parameter. The individual arcs can never be singular and are always convex. An inflection point of an A-spline must be placed at a junction point. © 1998 Elsevier Science B.V.
1. Introduction A-splines are piecewise algebraic curves. Sederberg proposed their use in geometric modeling (Sederberg, 1984); Bajaj coined the name. The A-splines discussed in this paper are curvature continuous. This property is equivalent to G2; see (Degen, 1988; Hohmeyer and Barsky, 1989; Garrity and Warren, 1991). We restrict here to cubic curves, since these have sufficient degrees of freedom to form G 2 splines, but few enough that we can give a geometrically intuitive interpretation for each of them. Two applications for splines are free-form geometric design and Hermite interpolation. We will discuss both of these applications for cubic A-splines. When using an A-spline for free-form design one begins with a control polygon as in Fig. 1; a sequence of * Corresponding author. E-mail: [email protected]. J E-mail: [email protected]. 0167-8396/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PH S01 6 7 - 8 3 9 6 ( 9 7 ) 0 0 0 3 1 - 9
262
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287 P
'p 0
Fig. 1. An A-spline.
triangles Ti such that two adjacent triangles have collinear edges. The spline interpolates joints Pzi, is tangent there to each segment Pzi-lPzi+l and has prescribed curvatures at each junction point. Each of these curvatures can be set numerically or relative to a default value. Each segment of the spline is convex, but zero curvature is allowed at any of its joints. In particular this means that inflection points may be placed where required and cannot occur otherwise. For example in Fig. 1, the curvature at P4 must be zero. Each segment can be controlled independently, keeping fixed the curvatures at its endpoints. The tangency and curvature interpolation conditions at the joints use six of the nine available degrees of freedom. The remaining three are managed as geometric handles for each individual segment. These are: an interior interpolation point B, the tangent line at B within limits imposed by the curvatures at the ends, and a tension parameter r. When cubic A-splines are used for high accuracy Hermite interpolation, the points P2i are chosen on a given curve, and then using the tangent lines at the Pzi the control polygon is formed around it. If inflection points are present they must be selected to be among the P2i. The curvatures at the P2i are computed and become additional input, as well as another point B in each triangle and its tangent. Using the results of (Degen, 1993, 1995) one can see that for each fixed value of the tension parameter r the A-spline we construct here approximates the original curve to order 8. We start with an overview of some of the previous work on G 2 splines, beginning with parametric splines. A construction of G 2 conic splines is possible, but they can not interpolate curvatures given at the joints and have only one (global) shape handle (Pratt, 1985; Farin, 1989). Pottman describes G 2 splines created from elements, each of which is composed of two conic arcs pieced together (Pottmann, 1991b). These curves have local control, but can never have zero curvature and have fewer shape handles than cubic splines. Cubic spline interpolation using parametric cubics as described by Farin produces curves that interpolate points but not tangents or curvatures (Farin, 1993). The choice of a knot sequence is required, whose effect on the curve is not well enough understood. A global system of linear equations must be solved.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
263
The direct G 2 cubic splines ('y-splines) described in (Farin, 1993) interpolate only the endpoints and tangents. An alternative construction for these curves is given in (Pottmann, 1991a). A local scheme for the use of polynomial cubic curves for high accuracy (order 6) Hermite interpolation was proposed by (de Boor, H011ig and Sabin, 1987). The method involves solving a system of nonlinear equations. Shirman and S6quin extend this construction with more stress on the application to free-form design (Shirman and S6quin, 1992). They describe G 2 splines created from elements, each of which is composed of two cubic arcs pieced together. Degen extends the construction to rational curves and achieves interpolation of order 8 (Degen, 1993, 1995). The most difficult part concerns the existence and reality of the curves. Sederberg (1984) proposed the use of algebraic splines and Warren (1991) expanded on his ideas. Li, Hoschek and Hartmann (1990) gave a construction that produces high order geometric continuity at each junction point but uses only one of the available degrees of freedom for shape control. Levin and Nadler (1995) extended these ideas to give methods for interpolating points that lie on the boundary of a convex region with an oval of an algebraic curve. Bajaj and Xu (1992, 1996) discuss tangent, curvature and higher order geometric continuity of A-splines. They use all of the available degrees of freedom (subject to one inequality), not for direct shape control, but for higher order approximation at the junction points and for least squares approximation of additional data points within each triangle. They also show how algebraic splines can be constructed in a sequence of parallelograms or rectangles (Bajaj and Xu, 1995). The authors studied the construction of cubic algebraic splines in (Paluszny and Patterson, 1993, 1994a), using all but one of the possible degrees of freedom for direct shape control. Both tangent continuous and curvature continuous splines were constructed, but the construction of curves that were not necessarily globally convex was hindered by not using all the available degrees of freedom. In this paper all of the available degrees of freedom are used for direct shape control. We find a larger family than Bajaj and Xu were able to use, and the restriction to globally convex curves is removed. In fact we characterize the maximal family of cubic curve segments that satisfy the interpolation, tangency, and curvature conditions at the given points. There have been two main objections raised against using implicit splines. First it was assumed to be impossible to control the shape of a curve in a geometrically intuitive way by adjusting the coefficients. This paper should lay to rest that objection; we show how to use all of the available degrees of freedom for geometrically meaningful shape control. Second it was asserted that rendering and robustness problems would never be solved. This objection is not dealt with in this paper, but some progress has been made; see the references cited in Section 9. Naturally, parametric curves are easier to render and methods for creating them extend to space curves and tensor product surfaces. By contrast, implicit curves provide easy control of convexity and displaying a family of solutions certainly solves the problem of existence and reality of solutions (among which will be any parametrizable curves that exist). The conclusion is to provide a hybrid system
264
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
that takes advantage of the strengths of both parametric and implicit curves (Marsilio et al., 1994). The organization of the paper is as follows. In the next section individual cubic segments are studied and the interpolation point B is introduced. Section 3 describes the parameter space for controlling the segment and introduces the slope parameter m. In Section 4 a cubic discriminant A0 is defined. The graph of A0 = 0 in the parameter space delineates the portion of the parameter space that corresponds to useful curves. This acceptable region is studied further in Section 5. Section 6 studies the tension parameter r. In Section 7 we explore the use of G 2 cubic A-splines for free form geometric modeling and some proposed defaults for the initial settings of the shape parameters. Section 8 discusses high order Hermite interpolation and Section 9 contains a brief discussion of graphing implicit curves. Notation: This paper uses three different types of coordinate systems. They are distinguished notationally as follows: Triples P[s, t, u] in square brackets refer to the homogeneous coordinates of the point P, triples in round brackets P ( s , t, u) to the barycentric coordinates of P, and pairs P ( x , y) in round brackets to Euclidean coordinates.
2. Behavior of a cubic segment
In this section we investigate the behavior of a cubic curve segment in an individual triangle To. The segment is defined by an implicit equation with ten coefficients; however Sederberg observed that the interpolation of endpoints and end tangencies in each triangle causes four of the ten coefficients to be zero (Sederberg, 1984). The equation for the ith segment can then be expressed in terms of the barycentric coordinate system defined by the ith triangle in the control polygon as F = asZu + bsu 2 - cst 2 - dt2u + e s t u - f t 3 = O.
The coefficients can be associated to trisection points of the triangle as in Fig. 2. The remaining six coefficients can be used for geometric control of the curve segment. One can think of the graph of F = 0 as the zero'th level curve of a functional triangular surface patch, but the multinomial coefficients have been absorbed into a, b , . . . , f to simplify notation. Some conditions on the signs of these coefficients are necessary. It is simple to check that a and b must have the same sign in order for the curve not to cross the segment POP2. We must make sure that the curve enters the triangle at/90. If we choose Euclidean coordinates so that the vertices are P0(0, 0), PI (mj, 0) and P2(x2, y2) (see Fig. 3) and compute the branch expansion of F = 0 in a neighborhood of P0, we obtain C Y2
---x Y=ax~
2
+.--.
From this it is clear that c and a must have the same sign for the curve to enter the triangle at P0. A similar result holds for coefficients b and d. Thus it is necessary to have a, b, c, d of the same sign to have any hope of finding an appropriate arc inside the triangle from P0 to P2, so we now impose this condition.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
265
P0/ s=l
1.1=1
Fig. 2. Association of coefficients to points of the triangle.
..... "'"'"'"'"'"'"'"'""~.0I"//~/xl, x
Po(O, o) Fig. 3. Branch expansion of F ----0 near 19o.
If in addition f is not zero and if the sign of f is the same as that of c and d the curve will not cross either PoP1 or P1P2. In this case the only boundary points of the triangle that the curve passes through are Po and P2. The branch expansion shows that Y=
c Y2 X 2
ax-~
is an osculating parabola to the curve at the origin. Its curvature at the origin is easily computed to be
~ = 2 cy2 a x2"
Hence in the notation of Fig. 4 we have the following result. Proposition 1. If F = 0 is an algebraic curve of degree 3 in which the coefficients of s 3 and 82t a r e O, its curvature at Po is given by
2 c hi a 9o~
4c_A a g o~
where A is the area of triangle PoP1P2. Similarly if the coefficients of u 3 and tu 2 are O, its curvature at P2 is given by dho =4 d A N,1 = 2 ~ if-'}ff-~.
266
M. Paluszny, R.R. Patterson / ComputerAided Geometric Design 15 (1998) 261-287
-,.
os S
Fig. 4. Euclidean measurements of the triangle. It turns out that having coefficients a, b, c, d, f all of the same sign guarantees that the curve will have an arc from P0 to P2 and no inflection point or singular point inside the triangle (although it can have a singular point somewhere outside the triangle). The following result is due to Bajaj and Xu (1992, 1996). Theorem 2. For a, b, e, d, f all o f the same sign, any straight line that passes through 191 and crosses the line segment PoP2 intersects the curve one and only one time (counting multiplicities) in the interior o f the triangle PoP~ P2. Furthermore the curve segment can have no inflection point in the interior o f the triangle. The following result generalizes part of the above theorem by removing the condition on the coefficient f .
Proposition 3. l f a cubic curve contains a nonsingular arc that enters the triangle at Po tangent to PoP1 and leaves at P2 tangent to PIP2, then there are no inflection points on this arc in the interior o f the triangle.
Proof. A real cubic curve has either one or three inflection points. In the case of three, they are collinear (Salmon, 1879). Since the cubic segment in the triangle has the same concavity at P0 and P2, if it has any inflection points in the triangle it must have two. The line through them either separates P0 and P2 or it does not. If it does separate/90 and P2 then after the curve crosses the line twice (see Fig. 5a) it is trapped on the P0 side and cannot get to P2 without another inflection. If the line does not separate P0 and P2 (see Fig. 5b), the center portion of the segment, which must be concave up, cannot be joined to P0 and P2 without crossing the line again. D We now begin the analysis of how the coefficients can be used for geometric control of the curve. We can use two of the remaining degrees of freedom to specify curvatures at the endpoints. Referring to Proposition 1, define #0 = c / a and #l = d/b. These can be used to control the curvatures at the joints. From the branch expansion, both #0 and/zl must be nonnegative. The values e = #0a and d ~1 b can be substituted into F. =
267
M. Paluszny, R.R. Patterson / ComputerAided Geometric Design 15 (1998) 261-287
(a)
(b)
P2
Fig. 5. Inflection points inside the triangle.
There are three remaining degrees of freedom to be specified. For the first of these, a point B(s0, to, u0) in the triangle is selected. Requiring the curve to interpolate B allows one to solve and determine a specific value e0 for the coefficient e: eo = - ( a s 2 u o + bsou~ - a # o s o t ~ -
b#lt2uo - f t 3)/sotouo.
With the additional substitution of e0 for e, the equation of F becomes asZu + bsu 2 - a#ost 2 - b#lt2u + eostu - f t 3 = O,
(1)
where #0, Pl, so, to and uo are now geometrically meaningful parameters. It remains now to provide geometric meaning for the remaining parameters a, f and b, and to do so in such a way as to obtain the largest possible family of useful curve arcs.
3. The [a, f, b]-parameter
plane
Many useful curves exist for which the sign of f does not agree with that of a and b. If f = 0 the curve passes through P1, and for f of opposite sign to that of a, b one gets a branch of the cubic coming into the triangle. This is not a problem so long as there is a segment from Po to P2 that is traced, say from Po to P2, and the two curve branches do not touch (the curve remains nonsingular inside the triangle). Bajaj and Xu do not allow these extra curves because they are using least squares methods to find the best fitting curve segment to data points inside the triangle--a significantly different problem. With #o, #1 and B chosen, the cubic depends homogeneously on the parameters a, f and b. One can think of the projective plane of [a, f, b] as the parameter plane for the family of cubics of the form (1) with curvatures determined by #0 at P0 and #l at P2 and interpolating the point B. The shaded region in Fig. 6 consists of those points A(a, f, b) for which ab > 0. By Theorem 2, at least every point in the triangle P a P f P b gives a cubic that has a convex segment connecting P0 and P2. Using B for the first shape handle, the other two arise from the following observation.
Proposition 4.
Cubic segments with the same #0, ~1 and B that have the same tangent line at B correspond to parameter points A that lie on a straight line in parameter space.
268
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 2 6 1 - 2 8 7
Fig. 6. Preliminary domain of A.
P,
, ~ ' ~ [ s 1, O, Ul] Fig. 7. Interpolating cubic arc. Proof. The equation of the tangent line to F at B is
Fs(B)s + Ft(B)t + F,~(B)u = 0. This line meets the PoP2 line at the point Q[sl, 0, Ul], where [Sl,0,?~l] = [ - Fu(B),O, Fs(B)]; see Fig. 7. The point Q lies outside the segment PoP2 because the cubic arc can have no inflection points inside the triangle. Such points Q are parametrized by m ~
81 q- U l 81 - - u 1
for-l~
--
1. Thus
81 -Jr Ul __ - -
81 -- '/3,1
Fu(B) - F,(B) f u ( B ) + F~(B)'
which simplifies to
asg [m(#oto2 + u~) - (#oto2 - u2)] + ft~ [rn(so + uo) - (so - uo)] + bug[m(p,tg + s~) + (#it 2 - sg)] = 0.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
269
Fig. 8. Preliminary strategy for controlling cubic segments. This is the equation of a straight line in the [a, f, b]-plane.
[]
The line in Proposition 4 is called an m-line. It is remarkable, but easy to check, that all of the m-lines pass through a common point
I[t30uo(souo- #,t~),souo(#o#,t 4 - s~u~),t30so(souo- #ot~)] • The point I exists except in the case in which/-to and #, happen to be the curvatures of a conic arc that solves the problem. Indeed, if the notation 7 = t~/souo is introduced, the coordinates of I are / [~o7(1 - #17), #o#,72 - 1, u~7(1 - #o7)] ;
(2)
and all three coordinates of I are zero in the case #0 = #1 = 1/7. In this case the cubic F is reducible:
F = (Tsu - t 2) (as + f T t + bu). The quadratic factor is a solution to the problem: it satisfies the interpolation conditions at P0 and P2, has the curvatures determined by #0 and #1, and interpolates B. The parameters a, f and b have no affect on the conic. Proposition 4 suggests a strategy for controlling the cubic segments: select a point B to be interpolated, select a tangent line at B by means of choosing an m-line in the pencil through I, and then for the third shape handle use the remaining parameter r to move the parameter point A along this line. The situation is much like polar coordinates; see Fig. 8. However, the limits on m and r can be extended beyond the PaPyPb triangle to a natural boundary that is given by a certain curve Ao = 0 in the parameter plane.
4. The cubic discriminant /t
There is a universal polynomial A of degree twelve whose variables are the ten coefficients of a general homogeneous cubic polynomial in three variables (Salmon, 1879).
270
M. Paluszny, R.R. Patterson / Computer Aided Geometric
Design 15 (1998) 261-287
It has the property that when d > 0 the cubic has one circuit (as a curve in projective space), when L < 0 the cubic has two circuits, and when n = 0 the cubic is singular or reducible. In the case under consideration, in which four of the coefficients are zero, the formula for n simplifies (Paluszny and Patterson, 1993). After the substitutions c = pea and d = ~1 b, A reduces further to A=a3b3(-e3f3+27abf4+e4f2po-36abef3po+8abe2f2& + 16a2b2f2& + e4f2p1 - 36abef3p1 - e5f,uopl +46abe2f2popul - 8abe3f&l
- 24a2b2f2&I
- 24a2b2f2p&
- 16a*b*ef&~
- abe4p&: + 64a*b*ef&T
+ 16a2b2f2& - 16a2b2efp&
+ 8abe2f2p: - 8abe3fp& - 8a2b2e2&T - 16a3b3&f
- 8a2b2e2p&: + 32a3b3&:
- 16a3b3&;)
= a3b3A0. The additional substitution e = ea creates the sixth degree homogeneous polynomial Ao(a, f, b), which is too long to reproduce here, but can easily be found with a computer algebra program. The graph of A0 = 0 in the [a, f, b]-plane will help to delineate a larger region in which the parameter points A(u, f, b) that determine useful cubic arcs may be selected. Proposition 5. When not both po and ~1 equal l/y, I never lies inside the triangle P,PfPb.
the point I is a singularity
of
do.
Proof. The coordinates of I satisfy A0 = 0, (no), = 0, (A,), = 0, and ( Ao)b = 0, so it is a singular point of multiplicity at least two. Using (2) one sees it is impossible for q all three coordinates of I to have the same sign, so I lies outside the triangle. The graph of A0 differs considerably depending on the relative sizes of po, pi and l/y. To simplify notation, introduce the variables va = ,LLOY - 1 and v1 = ply - 1 by making the substitutions PO = (uo + l)SO210& in Ao. In this notation I are
and
PI = (211+ ~)souo/$,
the formula for A0 becomes shorter again and the coordinates
1 [&or4, -(vov1 + uo + u,,&L;,
so&].
of
(3)
Because y > 0 and ~0, ~1 3 0, one has ~0, vi > - 1. This region in (~0, vl)-space is the meta-parameter space for studying Ao. The origin corresponds to the case in which I does not exist. In order to determine the type of the singularity I for A,, the Hessian of A0 must be examined as a quadratic form (Patterson, 1988a):
HeNAo)
=
(no),,(l)
(Ao)&)
(Ao)ab(I)
(AoLd
(Ao>.~J(~
(Ao)fb(r)
(Ao)ab(l)
(Ao),b(r)
(Ao)bb(l)
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
271
vI
(0,o)
Vo
(-1,-1)
Fig. 9. Graph of h and the meta-parameter space.
The following computation can be carried out using computer algebra. First each entry of Hess(Ao) is found to have the factor Sot 8,12 0 u o8-2 n , where 2 21 . h ( v o , Vl) = v 2 - 2 v o v , + v~ + 2 v 2 v , + 2 v o v 2 + VoV
Assuming that vo ¢ 0 and h ¢ O, the remaining matrix of Ao can then be reduced by a sequence of identical row and column operations to the matrix
(ooo) 0 0
h 0
0 0
.
Therefore the graph of h (see Fig. 9) separates the meta-parameter space into the region in which h > 0 and I is a crunode and the shaded region in which h < 0 and I is an acnode (isolated singularity). This result holds on the line Vl = 0 since no assumption was made about vl. It holds also on the line v0 = 0 since a corresponding result can be obtained by reducing the matrix using the (3, 3)-element of the Hessian as the first pivot instead of the (1, 1)-element. On the graph of h = 0, I is a third order singularity.
272
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
(b)
Fig. 10. Example graphs of Ao when vo, Vl > O and (a) vo = ~Ul
=
3, (b) vo = 3, vj = 4.2.
Here are some examples of the graph of Ao. In case v0, vl > 0, I is a crunode in the region ab > 0. In Fig. 10a, vo = Vl = 3. In this special case of vo = vl the sextic Ao factors into a quartic and a double linear factor whose graph is a double tangent. The point I is circled. In Fig. 10b, v0 = 3 and Vl = 4.2. In these examples the graph of F = 0 interpolates B ( 1 / 4 , 1/2, 1/4). In case v0, vl < 0, I also lies in the region ab > 0. In Fig. l l a , v0 = vl = - 0 . 5 . Again Ao reduces to a quartic and a double line. In Fig. 1 lb, vo = - 0 . 4 5 and vl = - 0 . 5 . In these examples (v0, vl) lies in the region in which I is an acnode and the graph of F = 0 interpolates the barycenter B ( 1 / 3 , 1/3, 1/3). Notice that the graph of A0 = 0 sometimes enters the PaPfPb-triangle. A curve F = 0 that corresponds to a point A on the graph of Ao = 0 in the triangle is singular, but its singular point does not occur
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
273
•
f
"--.-,~.~..,,,..~
(b)
®
Fig. 11. Example graphs of Ao when vo, vl < 0 and (a) vo = vj -- - 0 . 5 , vl = - - 0 . 5 .
(b) vo = - 0 . 4 5 ,
274
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
Fig. 12. Example graph of Ao when vj < 0 < vo.
inside the PoPl P2-triangle. These singular curves F = 0 must be the polynomial and rational solutions found by (de Boor, Hollig and Sabin, 1987; Degen, 1993, 1995). Fig. 12 shows one example in which v0 = 1 > 0, v~ = - 0 . 5 < 0, and B is the barycenter. Now I is a crunode that does not lie in the region ab > O.
5. The acceptable region Parameter points A that produce useful segments of the corresponding curves can be chosen to lie in the region ab > 0 below the PaPyPb triangle and above the graph of A0 = 0. (The previous sentence needs to be understood in the sense of projective geometry, in which the region below the triangle is connected to the region above the triangle through the line at infinity.) This region is called the acceptable region. By
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
275
1)l
4 1)0
Fig. 13. Areas in the meta-parameter space.
sampling parameter points in other regions one finds that the corresponding curves do not satisfy the requirements. The shape of the acceptable region depends on the point (v0, Vl ) in the meta-parameter space. This point is determined by the curvatures at the endpoints as reflected in #0 and #l and on the interpolated point B. Number the areas in meta-parameter space as in Fig. 13. For the case of area 1, in which v0, Vl > 0, the acceptable region is four sided with two curved edges and I at a vertex, as in Fig. 14a or 14b. The two curved edges in general do curve outward as shown, so all of the points in the region are visible from the point I. Most of them lie in the acceptable quadrilateral determined by I and the three vertices Pa, Pf and Pb. For the case of area 2a or 2b, in which v0, vl < 0, an arc of A0 = 0 crosses the lines a = 0 and b = 0, separating the point I from the triangle. In this case the acceptable region is three sided as in Fig. 15a or 15b. The curved edge in general does arch away from the triangle. Most of the acceptable points lie in the acceptable triangle CaPfSb (Fig. 15a) or SaPfCb (Fig. 15b). The remaining cases in which one of v0, vl is positive and the other negative are similar to the previous, except that the curve forming the lower edge of the three sided region passes through Pa or Pb, as in Fig. 16. Now imagine the m-line through I sweeping across the acceptable region. In many cases it is necessary to find an intersection of this line with the graph of A0 = 0. By Bezout's Theorem there can be as many as four intersections; finding the correct one would be possible but slow. For this reason it is proposed to restrict to the acceptable quadrilateral or the acceptable triangle. This compromise misses some acceptable curves but allows many more than the restriction to the original triangle PaPfPb. In order to find the limits on m in each case it is necessary to compute the value of m at each vertex of the acceptable region. This is easiest for area 1. In the v0, '/31 notation, the equation of an m-line is
a.o~o [m(.o~o * ~o + .o) -- ("o~o -- ~o + "o)] + St~ [m(~o + "o) + ~o -- ~o] + b~o~ [m(~o~, + ~o + ~o) + (~o~, + ~o - ~o)] = o.
(4)
276
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
Fig. 14a. In area 1, to cover the acceptable quadrilateral, m varies from the value m a ~-
80Vo -- UO + 80 sOVO + uo + 80
for the m-line through Pa to m b .~
'/z0v I ~ u 0
-- 80
UoVl -{- UO ÷ sO
for the m - l i n e t h r o u g h Pb. T h e n r v a r i e s so as to m o v e A ( a , f , b) f r o m I to o n e o f s e g m e n t s P a P I o r P I P b , d e p e n d i n g o n m (see Fig. 14a, b). T h e c u t - o f f v a l u e is
mf-
8o
-
the
u0
so + uo
which is the value for the m-line through Pf. This verifies the information in the first row of the table in Fig. 17.
M. Paluszn); R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
277
Fig. 14b. Acceptable region when vo, vl > 0 and (a) I lies below the triangle, (b) I lies above the triangle.
For areas 2, 3 and 4, since some of the vertices of the acceptable region are intersections of the graph of A0 = 0 with the lines a = 0 and b = 0, it is necessary first to find the coordinates of the points Sa, Sb, Ca, C6 and determine which ones are actually the vertices. When a = 0, A0 factors
-- f 74 Vl ( f Tto + buov l ) 3 ( f~/to + buo + b~tovl ) ( fa/tovo + buo v l + buovo v l ). The cubic factor corresponds to the point Tb[0, uovl,--')'to], where A0 = 0 has the line a = 0 as an inflectional tangent; see the graphs of Ao = 0. The next factor determines a point Sb[O, uo(1 + Vl),--Tto], and the last factor determines Cb [0, uovl (1 + vo),--Ttovo]. Similarly when b = 0, Ao factors
- f @ v o ( f T t o + asovo) 3 (aso + f~/to + asovo)(asovo + fTtovl + asovovl).
278
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
St)
$a Cb
!
I
Fig. 15a. The corresponding points are Ta[-Tto, sovo, 0], Sa[-Tto, so(1 + vo), 0]
and
Ca[-TtOVl, sovo(1 + vL), 0].
In order to distinguish the points SD, Cb in the graph of Ao one can use the cross ratio. Given four distinct points P/J1,0],/95[0, 1], U[u, 1], V[v, 1] on the projective line, their cross ratio is given by
cR(Ps, Pb, u, v) =
1
0
0
1
u
1
v
1
v
0 u
1 1
1 v
0 1
u
The points U and V are separated by Pf, Pb if and only if C R < O. If U and V are not separated by Pf, Pb then they occur in the order P f P b U V if v / u > 1 and in the order P f P b V U if v / u < 1.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
279
,I
Cb
Sa
Fig. 15b. Acceptable region when vo, vt < 0, I is an acnode, and (a) I lies below the triangle, (b) I lies above the triangle.
Applied to the problem at hand one has 1
0
CR(Ps P~,S~,C~)= uo(1 + v~) -Tto '
0 uo(1 + v l )
1 -~to
0
1
uovl(1 + vo) -Ttovo 1
0
u0vl(l + vo)
-~/tovo
(vo+ 1)v~ (v, +l)vo Thus the four points P f , Pb, Sb, C6 occur in this order if and only if 'vt < v0; in area 2a. Similarly the points Py, Pa, So, Ca Occur in this order if and only if v0 < vl; in area 2b. It follows that the vertices of the acceptable triangle are P f , Co and Sb in area 2a, and p f , S,~ and Cb in area 2b.
280
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
cb
Fig. 16. Acceptable region when v~ < 0 < vo.
Area
1
Values of v0, v~
Vertices
vo, vl > 0
P~,Ps,Pb, I
Domain of m m a to m b
2a
vl < v o < 0
C~,Ps, Sb
Tt~ca to ft2sb
2b
vo < vl < 0
Sa, Ps, Cb
r o s a to T/Zcb
3
vl < 0 < v o
P,~, P~ , Sb
m y to m s b
4
vo<0
S,~,Ps,Pb
rns~ t o r a f
Fig. 17. Limits of the acceptable region. Two more rows of the table in Fig, 17 are completed by finding the values of m for which the m-line passes through Sa, Sb, Ca, Cb. These are respectively rosa = - 1 , m s b ~- 1,
281
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
(b) Fig. 18. Region in which I lies in Case 3.
mca
=
(s0
-
uo)(vo
-
-
uovo
l
(SO + U0)(V0 -- Vl) + U0V0U1
and
mcb =
(so
-
u0)(
o -
-
sovo
l
(SO + U0)(V0 -- Vl) -- 80V0Vl
.
In the case of area 3, v0 > 0, vl < 0 as in Fig. 16, so the vertices of the acceptable triangle are Pa, Py and either Sb or Cb. Since both components of Cb have the same sign, Cb lies between P f and Pb. However the signs of the components of Sb are opposite, so 8b is the comer of the acceptable triangle. In Case 3 the only other thing needed to complete the table is to verify that the point I is located in the [a, f, b]-projective plane in such a way that the extreme vertices of the acceptable triangle as seen from I are always P I and Sb, never Pa. Note that the line t3of + (1 + Vl)Sou2b = 0 through Pa and Sb and the line b = 0 determine two regions in the projective plane, and that the triangle P ~ P I P b lies in one of them. It is enough to check that for v0 > 0 and vl < 0, I is contained in the same region as the triangle. Figs. 18a and 18b illustrate this for two possible locations of the point Sb. The computation for area 4 is analogous to that for area 3.
6. The tension parameter r Once #0, #1, B and m have been chosen, there is a family of curves F = 0 determined by the remaining parameter r. These curves have three points of contact at each of Po and P2 since they have the same tangent line and curvature, and two points of contact at B since they have the same tangent line. Because all cubics through eight fixed points pass through a ninth fixed point (Walker, 1950), it follows that all curves in this family share another point J. If J does not lie on the curve arc from Po to P2 inside the triangle, then for any two curves of the family, one lies entirely to one side of the other except at the endpoints and B. And in fact, most of the time r serves as a tension parameter, pulling the curve towards or away from the polygon determined by the sides u = 0, s = 0 and the tangent line at B.
Proposition 6. /f/z0, #1, B are fixed so that (vo, vl ) lies in area 1 or 2, and m is within the bounds delimited by the acceptable region, then the ninth point o f intersection o f any two cubics does not lie on the curve arc inside the triangle.
282
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
Proof. For (v0, V,) in area 1 or 2 let Q(ao, 0, bo) be the point where the m-line intersects the line P~Pb. Parametrize the m-line by (1 - r ) l + rQ. Choose two distinct values rl and r2 and let F1 and F2 be the corresponding cubics. The t-resultant of F1 and F2 is
]~(r 1
-
-
T2)3(VoVl -'~ V0 -Jr-U1)S3U3(S/Z0
-
-
2tSO)2(aoa2s + b0/32u),
where the factor k depends only on B, and a and/3 in the last linear factor LI are
a = aosov~ + bouov~ + bouovov~, /3 = aosov 2 + bouov~ + aosov2v,. The factors s 3 and u 3 correspond to P2 and P0, the factor (suo - uso) 2 corresponds to B, and the linear factor gives a line L, = 0 through Pl and J. Similarly the s-resultant of Fl and F2 has a linear factor L2 that represents a line through Po and J :
L2 = -aosouo(vovl + vo + vl ) a t + to(aosovo - bouovj )flu. Solving Ll and L2 simultaneously one finds the coordinates of J
d [bosouo(vovl + vo + vt)/32, to(bouov, - aosovo)a/3,-aosouo(vovl + vo + Vl)a2]. When m is such that Q lies between Pa and Pb, ao and b0 are positive, so the first and third coordinates of J have opposite signs and J cannot even lie in the region su > O, thus cannot lie inside the triangle PoP1P2. In particular when (v0, vl) is in area 1 and A is in the acceptable quadrilateral J cannot lie inside the triangle. However when m is varied so that Q no longer lies between P~ and Pb, J can enter the triangle. This situation is easier to examine using actual values of a0 and b0, which are determined from (4):
ao = uo(m(uovl + uo + so) + uov, + uo - so) bo = - s o ( m ( s o v o
and
+ uo + ~o) - ~ovo + ~o - so).
The coordinates of J become J[j,~, jr, ju], where
js = - s o ( m ( s o v o + uo + so) - sovo + uo - so)/3 2, j~ = t o ( m ( ~ 0 + ~0) + ~0 - so ) a / 3 ,
j~
=
- u o ( m ( u o v , + uo + so) + uov, + uo
-
s o ) a 2,
and where now a -- m ( ( s o + ~ o ) ( ~ o
-~l)
- so~o~l) - ((so - - o ) ( ~ o
- ~,) - ~o~o~,)
and /3 -- m ( ( s o + uo)(vo - v , ) + u o v o v , ) - ((so - uo)(vo - ~,) - u o ~ o v , ) .
Thus one sees that as m sweeps the acceptable region, js = 0 only when m = inca, which makes/3 ----0, or m = m~, which makes the other factor zero. For (vo, Vl) in area 2a, race is a bound of the domain of m, while in area 2b, mc~ is outside the domain of m. So for area 2, /3 can vanish only at one end of the domain of m. For m = m~
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
283
the coordinates of J are [0, t0vo(1 + v l ) , u o ( v o v i + vo + Vl)]. For (vo, vl) in area 2, jt and ju have the same sign, so J lies on the P1P2-side of the triangle. Thus J can enter the Pi P2-side of the triangle only on the other branch of the cubic. A similar argument holds for the PoPl-side of the triangle. So for (v0, vl) in area 2, J cannot lie on the arc from P0 to P2. The only case that remains to be proved is when (v0,vl) lies in area 1 and A is in the acceptable region but outside the acceptable quadrilateral. This proof is omitted because of the decision not to use these curves anyway. [] So r acts as a tension parameter whenever both p0,,ttl > 1/7 ((U0,'UI) lies in area 1) or #0, #1 < 1/7 ((v0, vl) lies in area 2). For some choices of (v0, vl) in areas 3 and 4, J does lie on the curve arc from P0 to P2 inside the triangle and can even be B. In practice in these cases it is observed that motion of A along such an m-line is essentially doing very little to change the graph of F = 0 anyway. Here are some illustrations of the use of the tension parameter r. For v0, vl > 0 (area 1), when A approaches one of the bottom curved boundaries of the acceptable region, the curve develops a comer as in Fig. 19a. In the extreme case of A actually on the bottom boundary the curve has a crunode in triangle PoP1P2; when A = I the crunode is at B. For vo, Vl < 0 (area 2), when A approaches the bottom curved boundary of the acceptable region the curve fits tightly against the tangent line at B as in Fig. 19b.
7. G 2 cubic A-splines in free-form m o d e l i n g
In this section we discuss some defaults that can be used so that as soon as a designer has constructed a control polygon as in Fig. 1, a curve will appear. First consider curvature. Curvature is an invariant of Euclidean geometry and introducing it is somewhat unfortunate since all the other constructions with F are being done in affine geometry. Without knowing the scale of a figure it is impossible to appreciate the significance of a number as curvature or radius of curvature at a point. One consequence of having defaults for the curvatures is that a slider can then be used to control the curvature at a junction point. In order to construct default curvatures at the junction points, consider two adjacent triangles T i - i and Ti and the graph of su - t 2 = 0 in each. In the notation of Fig. 4, the curvature of the segment in Ti-l at P2i is 2hzi_z/g2i_l while that of the curve in Ti is 2h2i+1/92i. Ordinarily these curvatures will be different. However if we define )ki--
h2i-2g~i 2 h2i+lg2i_ 1
and use V/~iisu - t 2 = 0
in Ti-1
and
- - 1s u
_ t2 = 0
in Ti,
284
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
(a)
(b)
Fig. 19. Example cubic arcs for A near the boundary of the acceptable region. these conics both have the same curvature at P2~, namely the geometric mean of the original two curvatures ~i = 2 ~ / h 2 i - 2 h 2 i + l g2i-lg2i
Although we cannot generally create a G 2 c o n i c spline in this way, we can use ~i as the default curvature at P2i for a cubic spline. If the spline is open, at the endpoints we can use the curvature of s u - t 2 = 0; hi
~O=2g-~
2
and
~'~=
h2n-2
g22n_l"
The next defaults are the interior points B~ that are interpolated; they can be chosen to be the barycenters. Choosing B~ also has the effect of selecting an acceptable region
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
285
in parameter space for control of the ith segment of the spline. The limits on the slope parameter mi at B~ are determined by the extreme vertices of the acceptable region as seen from I, whose position depends on Bi and the curvatures at the endpoints of the ith segment; see Fig. 17. The parameter mi can also be controlled by a slider, with the default being an intermediate value. The tension ri can then be controlled by yet another slider that moves the point A on the segment of the m-line in the projective plane. The limits on ri are determined by the value of mi and the line segments that bound the acceptable quadrilateral or triangle, and like m~ it can default to an intermediate value.
8. G2 cubic A-splines for Hermite approximation Suppose a smooth curve C is given, either parametrized or at least in such a way that it is possible to select points on it and compute tangents and curvatures. Choose points P2i on C, including any inflection points. The tangents at the P2i determine a sequence of triangles. Choose another point Bi on C in each triangle. We have seen how to construct cubic segments from P2i through B to P2i+2 that have the same tangent and curvature at each end, the same tangent at Bi and one remaining degree of freedom r. Let Mi be the midpoint of P2i and P2i+2. Because both C and the cubic segment (for any fixed r) are convex, rays from M~ provide a transverse direction along which distance can be measured between the two curves. Degen has shown that the maximum of this distance goes to zero as the eighth power of the distance h from PEi to P2i+2 as the partition is refined (Degen, 1993, 1995). (Actually he measures distance h in a parameter space, but this is equivalent.) This approximation order could theoretically be increased to nine if a parameter could be found to replace r and increase the order of contact at B or if the curve were to interpolate two interior points instead of just one. If a curve segment comes from a conic it will be reproduced exactly; cubic A-splines have the property of quadratic precision. Because of the endpoint and tangency conditions on the conic, in its equation the coefficients of s 2, '/.Z2, 8 t and tu will be zero. We can then normalize the equation to the form 7su - t 2 = 0. Because B lies on the curve, = t~/(souo). The curvatures of this conic at Po and P2 are given by n0 -
2 hi
"y go
and
/~1 :
2 ho
---
gF
So /~0 :
~1 :
l/"y.
This is the case in which I does not exist and the cubic is reducible; see Section 3. The conic is one of the factors and is the arc that is traced from P0.
9. Graphing implicit curves There are a number of methods available for graphing implicit curves: (Bajaj et al., 1988; Chandler, 1988; Patterson, 1988b; Sederberg et al., 1989; Dobkin et al.,
286
M. Paluszn3, R.R. Patterson/ Computer Aided Geometric Design 15 (1998) 261-287
1990; Bajaj and Xu, 1994; Marsilio et al., 1994; Paluszny and Patterson, 1994b; Taubin, 1994). The simpler of the two methods described by Taubin was used to graph h and A0 in Figs. 9, 10-12, 14-16. Cubic A-spline segments can be traced, since their behavior inside each triangle is so predictable. The fastest and most robust method that we have found to graph them uses a ray that emanates from the midpoint of the PoP2 segment and rotates in steps from P0 to P2. The ray is intersected with the tangent line to the curve at the previously computed point, and then a Newton correction is made to find the next point on the curve. This method was used for Fig. 19.
Acknowledgments The authors wish to thank W.L.E Degen for reading a draft of this paper and providing many helpful suggestions. This research was partially supported by the NSF/Conicit U.S.-Latin America Cooperative Science Program, grant INT-9304208.
References Bajaj, C., Hoffmann, C., Hopcroft, J. and Lynch, R. (1988), Tracing surface intersections, Computer Aided Geometric Design 5, 285-307. Bajaj, C. and Xu, G. (1992), A-splines: local interpolation and approximation using Ck-continuous piecewise real algebraic curves, Computer Science Technical Report, CAPO-92-44, Purdue University. Bajaj, C. and Xu, G. (1994), Rational spline approximations of real algebraic curves and surfaces, in: Dikshit, H.P. and Micchelli, C.A., eds., Advances in Computational Mathematics, C, World Scientific Publishing Co, Approximations and Decomposition Series, in press. Bajaj, C. and Xu, G. (1995), Regular algebraic curve segments (III) - Applications in data fitting, manuscript, Purdue University. Bajaj, C. and Xu, G. (1996), Data fitting with cubic A-splines, in: Gigante, M. and Kunii, T., eds., Proceedings of Computer Graphics International, CGI'94, Melbourne, Australia, World Scientific Publishing Co. Chandler, R. (1988), A tracking algorithm for implicitly defined curves, IEEE Computer Graphics and Applications 8, 83-89. de Boor, C., HOllig, K. and Sabin, M. (1987), High accuracy geometric Hermite interpolation, Computer Aided Geometric Design 4, 269-278. Degen, W. (1988), Some remarks on Brzier curves, Computer Aided Geometric Design 5, 259-268. Degen, W. (1993), High accurate rational approximation of parametric curves, Computer Aided Geometric Design 10, 293-313. Degen, W. (1995), High accuracy approximation of parametric curves, in: Dahlen, M., Lyche, T. and Schumaker, L.L., eds., Mathematical Methods in CAGD II1, to appear. Dobkin, D., Levy, S., Thurston, W. and Wilks, A. (1990), Contour tracing by piecewise linear approximations, ACM Transactions on Graphics 9, 389--423. Farin, G. (1989), Curvature continuity and offsets for piecewise conics, ACM Transactions on Graphics 8, 89-99. Farin, G. (1993) Curves and Surfaces for Computer Aided Geometric Design, Academic Press, New York.
M. Paluszny, R.R. Patterson / Computer Aided Geometric Design 15 (1998) 261-287
287
Foley, T., Goodman, T. and Unsworth, K. (1994), An algorithm for shape preserving parametric interpolating curves with G 2 continuity, in: Laurent, P.-J., Le Mrhautr, A. and Schumaker, L.L., eds., Curves and Surfaces in Geometric Design, A K Peters, Wellesley, MA, 249-259. Garrity, T. and Warren, J., (1991), Geometric continuity, Computer Aided Geometric Design 8, 51-65. Hohmeyer, M. and Barsky, B. (1989), Rational continuity: parametric, geometric and Frenet frame continuity of rational curves, ACM Transactions on Graphics 10, 335-359. Levin, D. and Nadler, E. (1995), Convexity preserving interpolation by algebraic curves and surfaces, Numerical Algorithms 9, 113-139. Li, J., Hoschek, J. and Hartmann, E. (1990), G n ~-functional splines for interpolation and approximation of curves, surfaces and solids, Computer Aided Geometric Design 7, 209-220. Marsilio, R, Paluszny, M., Patterson, R. and Rodriguez, O. (1994), Ambiente h~rido de disefio gr~ifico: splines algebraicos y paramrtricos, in: Memorias XX Conferencia Latinoamericana de Informdtica, Mexico City, Mexico, 147-156. Paluszny, M. and Patterson, R. (1993), A family of tangent continuous cubic algebraic splines, ACM Transactions on Graphics 12, 209-232. Paluszny, M. and Patterson, R. (1994a), A family of curvature continuous cubic algebraic splines, Technical Report #94-8, Department of Mathematical Sciences, Indiana University - Purdue University Indianapolis. Paluszny, M. and Patterson, R. (1994b), G2-continuous cubic algebraic splines and their efficient display, in: Laurent, P.-J., Le Mrhautr, A. and Schumaker, L.L., eds., Curves and Surfaces in Geometric Design, A K Peters, Wellesley, MA, 353-359. Patterson, R. (1988a), Parametric cubics as algebraic curves, Computer Aided Geometric Design 5, 139-159, Patterson, R. (1988b), Parametrizing and graphing nonsingular cubic curves, Computer-Aided Design 20, 615~i23. Pottmann, H. (1991a), A projective algorithm for curvature continuous rational splines, in: Farin, G., ed., NURBS for Curve and Surface Design, SIAM, New York. Pottmann, H. (1991b), Locally controllable conic splines with curvature continuity, ACM Transactions on Graphics 10, 366-377. Pratt, V. (1985), Techniques for conic splines, in: Proceedings of SIGGRAPH 85, ACM, New York, 151-159. Salmon, G. (1879), Higher Plane Curves, Hodges, Foster and Figgis, Dublin. Sederberg, T. (1984), Planar piecewise algebraic curves, Computer Aided Geometric Design 1, 241-255. Sederberg, T., Zhao, J. and Zundel, A. (1989), Approximate parametrization of algebraic curves, in: Strasser, W. and Seidel, H.-P., eds., Theory and Practice of Geometric Modeling, Springer, 33-54. Shirman, L.A. and S6quin, C.H. (1992), Procedural interpolation with curvature-continuous cubic splines, Computer-Aided Design 24, 278-286. Taubin, G. (1994), Rasterizing algebraic curves and surfaces, IEEE Computer Graphics and Applications 14, 14-23. Walker, R. (1950), Algebraic Curves, Dover, New York. Warren, J. (1991), Several notions of geometric continuity for implicit plane curves, Monograffas de la Academia de Ciencias de Zaragoza 2, 121-129.