- Email: [email protected]
Preserving convexity using piecewise cubic interpolation
Comput.&GraphicsVol. 15. No. I, pp. 15-23, 1991
0097-8493/91 $3.00+ ,00 ~' 1991 Pergamon Press plc
Printed in Great Britain.
Technical Section
P R E S E R V I N G C O N V E X I T Y U S I N G PIECEWISE CUBIC I N T E R P O L A T I O N K. W. BRODLIE and S. BUTT School of Computer Studies, University of Leeds, Leeds, UK Abstract--This paper considers the problem of drawing a smooth cubic through a set of data points (x, i = 1, 2. . . . . N, where the y-values are dependent on the x-values--a common problem in scientific visualisation. A typical approach is to estimate the slope of the curve at each data point and construct a pieeewise cubic interpolant which is then easy to plot. An additional requirement is often to preserve the inherent shape of the data. While monotonicity can be preserved by suitable slope selection, it is known that the same cannot be achieved in general for convexity. However, this paper shows that by allowing (if necessary) two cubic pieces in some data intervals rather than one, convexity can always be preserved. Two methods are presented and illustrated with example data.
Yr),
I. I N T R O D U C T I O N
A fundamental problem in computer graphics is the drawing of a smooth curve through a set of data points (xr, Yr), i = 1, 2 . . . . . N. In many applications, particularly in scientific visualisation, the y-values are dependent on the x-values, and it is natural to seek an interpolant of the form:
y =y(x) which matches the given data. An attractive class of interpolating functions are piecewise cubic polynomials which are easy to construct and easy to plot. A typical procedure is to estimate the slope dr of the curve at each data point by some method; this then uniquely defines the cubic piece between each pair of data points. Many slope estimation methods have been suggested--a simple one is the osculatory, or cubic Bessel[l], method in which the slope at a data point is taken as the slope of a quadratic passing through the data point and its neighbour on either side. A full review of different methods is given in Brodlie[5, 6]. The methods all have first derivative continuity at the data points by construction. Attention has focussed recently on methods that preserve some shape inherent in the data. For example, it may be known that the function underlying the data is monotonic, and it is desirable to impose monotonicity as a constraint on the interpolant. The monotonicity case was explored by Fritsch and Carlson [ 8 ] who proved the following result: suppose the underlying function is monotonically increasing (i.e., Y l <- Y 2 " • • -< YN), then the cubic interpolant in each interval [xr, x~+~] is monotonically increasing if the slopes at the end-points satisfy
the above inequalities. This means that given any monotonic data it is possible to construct a monotonic piecewise cubic interpolant. The Brodlie formula is used in the N A G Graphical Library [ 11 ] curve drawing routines, and has proved visually acceptable. Another shape property of importance is convexity which is the subject of this paper. Here the situation is less satisfactory. In contrast to the monotonic case, there exist convex data sets (i.e., data from an underlying convex function) for which it is not possible to construct a convex piecewise cubic interpolant. This is discussed further in Section 2, together with various solutions to the problem generally involving a different class of interpolant to handle the impossible cases. In Section 3 we show that convexity can be achieved with piecewise cubics; the impossible cases are handled by using two cubic pieces in an interval between data points. In Sections 4 and 5 two specific methods are derived: one ensures second derivative continuity at the join between the two cubic pieces; the other ensures a smooth transition from the case where a single cubic suffices to the case where two pieces are required. Section 6 shows the methods applied to a number of test datasets that have appeared in the literature. 2. PRESERVING CONVEXITY IN INTERPOLATION Suppose the data points are realisations of some underlying convex function. Thus Ar+l >-Ai
i = 1,2 . . . . .
N-2.
(2.1)
Suppose we assign slopes dr to the data points. In the interval [xr, xr+l ], the cubic polynomial interpolant is convex [ 12 ] if and only if
di < 3Ai
2d; + d,+~ _< 3A; _< d; + 2di+~.
(2.2)
dr+l -< 3A,,
Handscomb (Brodlie[6]) points out that this implies, for five consecutive data points,
where Ai = (Y;+I - Y i ) l ( x r + l - xr). Butland[3] and Brodlie[6] give simple formulas for dr in terms of(xr, Yr) and its two neighbours that satisfy
Ai+ 2 --
15
Ai- I ~
1 . 5 ( A i + I -- A i )
(2.3)
K. W. BRODLIEand S. BUTT
16
which is a condition on the data--thus there will exist convex data for which there is no piecewise cubic interpolant. Therefore, there is no incentive to seek a "magic" formula that will generate slopes di always satisfying the convexity condition (2.2), as was possible in the monotonicity case. Various strategies for preserving convexity have been suggested in the literature. Schumaker [ 13 ] shows how quadratic spline interpolation can preserve convexity: slopes di are estimated by any method one chooses, with the restriction that d~+l > dj. Then in each interval [x~, xi+t ], a quadratic spline interpolant is calculated, having at most one internal knot. Schumaker gives conditions on the position of that knot that ensure that the slope d* at the knot satisfies:
di+j >- d* > di. The resulting quadratic spline is necessarily convex. The drawback of the method is the fact that it has lower order than cubic interpolation methods. Another approach is to use piecewise rational functions. Gregory and Delbourgo[10] and Brankin and Gladwell [ 4 ] both suggest the use of piecewise rational cubics. The drawback of this approach is the additional complexity of plotting a piecewise rational cubic. Indeed Brankin and Gladwell only recommend the rational cubic when the simple cubic fails to be convex. Still another approach is to use parametric curves: Goodman and Unsworth [9] describe an interpolation scheme using parametric cubic polynomials that can generate single-valued curves that preserve convexity. However, a disadvantage of using parametric curves for the single-valued case is that it is relatively difficult to evaluate y for a given value ofx. 3. PRESERVING CONVEXITY BY USING T W O CUBIC PIECES
xl = 0 ,
y~ = 0 ,
x2=l,
y2 = 1,
with estimated slope at xt as dj = 0.1. Then, from the convexity condition (2.2), we see that the cubic Hermite interpolant is convex provided 1.45 < d2 < 2.8. For values of d2 outside this range, the cubic Hermite is nonconvex, and Fig. 1 illustrates the situation for d2 =5.
Theorem 1 Let (x~, yj ), (x2, Y2) be two data points, and let d~, d2 be the (estimated) slopes at these points. Suppose the data is convex i.e., dt < a < d 2 but the associated cubic polynomial interpolant is nonconvex i.e., d2 ~ [ t ( 3 A - d l ) , 3A -- 2d,] where A = Y2 - Yl/X2 Xl. Then we can insert a knot at x* E [xl, x2] and associate with that knot a value y* and slope d*, such that the piecewise cubic interpolant fitted over the two intervals [xl, x* ], and [x*, x2] is convex. Moreover the range of knot positions and associated values and slopes is given by: -
-
Suppose the data are convex, and we use a slope estimation method which ensures 1
di+l>di
i = 1,2 . . . . . N - 1 . 0.8
(Note that every common method will do this.) If condition (2.2) is satisfied, then the cubic interpolant in [x~, x~+l] will preserve convexity. So we need only be concerned with the case where (2.2) is not satisfied. The resulting cubic in this case is nonconvex, and we need some way of incorporating an extra degree of flexibility. We would prefer to do this while still using piecewise cubics. Dougherty, Edelman, and Hyman[7] suggest the relaxing of first derivative continuity at the data points, but this could yield a curve that is not visually acceptable. Our suggestion is to create the flexibility by the insertion of an extra knot between [xi, x~+~]. At this stage we simplify the notation by considering only the interval [x~, x2]. It will be useful also to have a specific example to consider in parallel with the theoretical development. We consider the case:
0.6
0.4
0.2
-012
O4
06 /
O8
-0.2
Fig. 1. Nonconvex cubic hermite interpolant.
t
Preserving convexity A, ~ [ m a x { d ~ , A -
](d2-
17
and so we extend ( 3.1 ) to
d,)},
m i n { A , d 2 - ~ ( d 2 - d~)}]
d~ < A, < min { A, d2 - ~ ( d2 - d, ) } .
(3.11)
A2 E [ m a x { A , A, + ~ ( d 2 - d , ) } , min {d2, A~ + ](d2 - d , ) } l d* E [max{ ½(3A~ - d,), 3 A 2 - 2d2},
( N o t e d~ < d2 - ~(d2 - d j ) , so (3.11) is c o n s i s t e n t . )
Next consider inequality (3.9). This can be written as:
m i n { 3 A , - 2d~, ½(3A2 - d2)}] where A, = (y* - y l ) / ( x * - x l ) , A2 = (Y2 - y * ) / ( x 2 -- X*).
A2 < A, + ](d2 - d~) and so we extend (3.10) to
Proof First, A, and A2 must satisfy: m a x { A , A, + ~ ( d 2 - d,)} < A2 d, < A, < A
(3.1)
A < A2 < d2
(3.2)
< min{d2, A, + ] ( d 2 -
dr)}.
(3.12)
This requires for consistency that Next d* must satisfy: A-
(3.3)
½(3AI -- d l ) < d* < 3A, - 2d,
(3.4)
½(3A2 - d*) -< d2 -< 3A2 - 2d*.
(3.5)
or
Al>A-~(d2-dl)
and so we extend ( 3.11 ) to max{ dl, A - ~(d2 - d l ) } , -< AI <-min{A, dz-~(d2-d,)} If inequalities ( 3 . 1 ) - ( 3.5 ) are satisfied, then the piecewise cubic interpolant is convex. Inequality (3.5) can be written as: 3A2 - 2d2 -< d* < ½(3A2 - d2)
(3.6)
and inequalities (3.3), (3.4), and (3.6) can be combined into the single inequality
which is easily shown to be consistent. Hence there exists a A~ satisfying (3.13), a A2 satisfying ( 3.12 ) and d* satisfying (3.7), and the theorem is proved. Note that the position of the knot (x*, y*) will be given by
(3.7)
using the fact that AI < ½(3Al -- d,) and A2 > ½(3Az - d2) which is true ifAl and A2 satisfy (3.1) and (3.2) respectively. Again if (3.1) and (3.2) are true, there will be a value of d* satisfying (3.7) provided that: ½(3A, - d,) < ½(3A2 - dE)
(3.8)
3A2 - 2d2 < 3Al - 2dl.
(3.9)
1
x*
- {(x2A2 - x l A l ) - - (Y2 -- Yl)} A2 - AI
y*
- { AiA2(x2 - Xl) - (A, y2 - A2yl)}. A2 - A3 (3.15)
max{ ½(3A, - d,), 3A2 - 2d2} ~ d* -< min{ 3A, - 2d,, ½(3A2 - d2)}
(3.13)
(3.14)
1
We can illustrate this theorem in terms of our worked example. Consider a slope d2 = 5.0 which causes the convexity test to fail. The theorem states that a knot can be inserted with associated slope d*, such that
+ )min(10 ,O;)]
Consider inequality (3.8) first. This can be written as i.e.,
A~ E5 [0.1, 1.01,
A2 ~ A, + ](d2 - d,) A2 ~ [ m a x ( 1.0, A, + - ~ ) , min(5.0, A, + ? ) ]
and we can extend (3.2) to m a x { A , A, + ](dz - d,)} -< Az ~ d2.
(3.10)
This requires for consistency that A,+ ](d2-d,)
or
A, < d 2 - ] ( d 2 - d,)
i.e.,
A2E[A,+-~,A,+?],
d* E [max{ ½(3At - 0.1), 3A2 - 10.0}, min{3A~ - 0.2, ½(3A2 - 5.0)}].
K. W. BRODLIEand S. BUTt
18
In Fig. 2, we choose A, = 0.5, A2 = 2.5, and d* = 1.0 and draw the resulting piecewise cubic which is convex. A knot is inserted as shown at (0.75, 0.375). Theorem 1 shows there is a range of possible knots and associated slopes. The following corollary shows that a simple choice of A,, A2 always lies within the range given in the theorem. Corollary l: Let A, = a A + (1 - a ) d i ,
A2 = a A + (1 - a)d2.
we require aA+(l
-a)d,
> A-
2(d2-dt)
i.e.,
2(d2-dl)>(1-a)(A-d,)
which is true for a E [ }, 2 ] since d2 > A. To show A , <--d 2 - -
]1( d 2 - d , )
we require
If a E [~, 2 ], then A1, A2 lie within the range given in Theorem 1 for all dr, d2, A, and the corresponding range for d* becomes:
aA + (1 -- a ) d , < d2 - ](d2 - d,) i.e.,
](d2-d,)>a(A-d,)
which is true for a E [0, ] ] . Thus
d* ~ [ m a x { d , + 3a(A - dj ), d2 - 3a( d2 - A)},
A, e [max{d,, A - ](d2- d,)},
min {d, + 3a(A - d,), d2 - 3a(d2 - A)}].
min{ A, d2 - ](d2 -
d,)}l
Proof." We have to show that
f o r a E [I, 21. Next we have to show that A, E [ m a x { d , , A -
2(d2 - d,)}, A2 E [max{A, A, + l ( d 2 -
min{ A, d2 - l ( d 2 - d , ) } ] .
d,)}, min {d2, Al](d2 - dj)}]
The results:
The results A, >_ d~
and
A, -~ A
A2> A
are trivial. To show A,>A-
and
A2
are trivial. To show A2 > A, + l ( d 2 - d , )
2(d2-d,) we require aA+(1
-- o~)d2 > aA + (1 - a ) d , i.e..
0.8
+ ~(d2-d~)
(1 - a ) ( d 2 - d , ) ~
~(d2-dt)
which is true for a E [ ~, 1]. To show A 2 <_ A, + 2(d2 - d , )
0.6
we require 0.4
otA + (1 - a)d2 -< aA + (1 -- a ) d , + ~(d2 - dr) i.e.,
(1-a)(d2-d,)_<2(d2-d,)
0.2
which is true for a E [0, 2]. Thus °
-0.2
0.2
0.4
0.6
0.8
1
A2 ~ [max{A, A, + -~(d2-
d,)}, min{d2, A,~(d2 - d,)}]
-0.2
Fig. 2. Convexity preserved by inserting a knot.
for a E [I, 21.
Preserving convexity The range for d* follows by simple substitution. Note: The extreme choices a = ] and a = 2 lead to particularly simple ranges for d * , viz..
a = ~ --,- d* ~ [ m a x { ' i ( d l + A), A}, min{A, ½(A+d2)}] a = ]--~d*@[max{A,
i.e.,
d* = A
2A-d2},
19
C o r o l l a r y 2: Let Al = a A + ( l - a ) d l ,
A2 = a A + ( 1
- a)d2
as in Corollary 1, with cr E [ -~, -] ]. Then the required value of d* for second derivative continuity is d* = [(d, + d2) - A] + 3a[A -- ½(d, + d2)]
min{2A--d1,A}]
i.e.,
d * = A.
4. A M E T H O D ENSURING SECOND DERIVATIVE CONTINUITY AT KNOT
To obtain a particular algorithm that will maintain convexity, we need to make a specific choice of A~, A2, and d* satisfying Theorem 1. It is convenient to work with the simpler formulas of Corollary 1. An attractive property is that the two cubic pieces be joined smoothly at the knot x * i.e., we seek second derivative continuity at x*, in which case the curve is a clamped spline. We have the following theorem.
We now explore for which values of a the value of d* given in Corollary 2 lies with the valid range of d* as given by Corollary 1. We have the following results. Theorem
3
Let Aj, A2 be as in Corollary 1. Then the value of d* giving second derivative continuity as identified in Corollary 2 lies within the valid range for all possible dl, d2 and A if and only if a = -~. P r o t f : Suppose a = ~. Then d* = A and it has already been noted after Corollary 1 that this is valid. Suppose a = ~ - e, for some { > 0. Then
Theorem 2
Consider a piecewise cubic c ( x ) constructed over the interval [x~, x2] so that
d* = a -
3~{a-
½(d, + d2)}.
Corollary 1 requires, among other things, c ( x l ) = )'j,
c(x2) = Y2,
c'(xl)
= d~,
c ' ( x 2 ) = d2
and with a knot at x* where x~ < x* < x2. Then the piecewise cubic has a continuous second derivative at x* if and only if
d*>_d2-3a(d2-A)=A-(1-30(d2-A). This requires: 3~(d2 - d l ) < d2 - A.
d ~
--
2 (x: - x, )
{3(hA 2 + k A , ) - ( h d 2 + k d j ) }
where h = x* - xl, k = x2 - x*. P r o o f Let c ~ ( x ) be the cubic piece over [x,, x * ] and c2(x) the cubic piece over I x * , x2]. Then it is easy to show (see Ahlberg, Nilson, and Walsh[2]) that 2dj 4d* c'~{x* ) = T + T 2(x
) -
4d*
60'2- y*)
k
k
k2
Method A
In any interval where the cubic polynomial is nonconvex, insert an extra knot with position determined by
6(y* - Yl) h~
2d2
Given any ~ > 0, we can always choose d~ < A < d2 such that this inequality fails. Hence the theorem is proved. This leads to the following method.
A, = 2A + ~d,
(4.1)
A2 = -~A + ~d2
(4.2)
Thus for c'f(x*) = c~(x*), we require
t~- + 4d*
+
-~ k _ 6(y* - yl) 4 6(y2 - y * ) h2
i.e., d * - - -
1
2(h + k)
and assign a slope d* = A at the knot. The resulting piecewise cubic polynomial is convex and is a cubic spline. Substituting (4.1) and (4.2) into (3.14) and (3.15) give the knot position as
k2
1
{ 3(hA2 + kay) - ( h d 2 + k d , ) } .
Theorem 2 tells us the following: once we have selected a knot x* and associated y*, we must choose a specific value of d* to obtain second derivative continuity.
X*
- {(x2d2 - x , d , ) d2 - dl
- (y2 - yl)}
y* = - [{d, d2 + 2A(d, + d2 - a ) } ( x 2 - x,) d2 - dl - 3 ( d , y2 - d2y,)].
20
K.W.
BRODLIEandS.
1
BUTT
1
0.8
0.8
0.6
0.6
0.4
0.4 0.2
0.2 i
0.2
-0.2
0.4
0.6
0.8
1
0.2
-0.i
Fig. 3. Method A--preserving 2nd order continuity at knot.
Figure 3 shows the effect of applying this method to our worked example. A knot is inserted as shown at (x*, y * ) = (0.82, 0.57) with slope d* = 1.0. 5. A METHOD INSENSITIVE TO CHANGES IN DATA Method A suffers from two difficulties. Firstly, it is a desirable property of any curve drawing method that the generated curve should not be sensitive to small changes in the data. Unfortunately this is not the case with Method A as the data moves from allowing a single convex cubic to requiring the addition of an internal knot. This is illustrated in Fig. 4, where the
0.4
0.6
0.8
1
Fig. 5. Always inserting a knot.
curves for the cases dE = 1.451 and d2 = 1.449 are shown alongside each other. This difficulty can be overcome by always inserting a knot, even when a single cubic will suffice. That is, we calculate Al and A2 by (4.1) and (4.2), and take d* = A, and use a pair of cubic pieces whatever the data. Since the expressions for At, A2 and d* are continuous with respect to the data, the resulting interpolant will likewise be continuous with respect to the data. This is illustrated in Fig. 5, where the curves for the cases d2 = 1.451 and d2 = 1.449 are seen to lie almost on top of one another. Contrast this with Fig. 4. However, there is a second difficulty, and this is the visual quality of the curve. The choice a = ] means that Aj and A2 are both relatively close to A, and so 1
0.8
0.8
0.6
/
0.4
,,," 06
0.4
0.2
o.~
o.',
o18
o'.8
;
0.2
C l a m p e d cubic spline Single cubic Fig. 4. Sensitivity to change of data.
OI
0.2
0.4
0.6
0.8
Fig. 6. Insensitive to change of data.
1
Preserving convexity the inserted knot is close to the chord joining the two data points and a "flat" curve results. In this section, therefore, we discuss an alternative method which is not sensitive to small changes in the data, and which aims to give a more visually pleasing curve. Second derivative continuity at the knot is sacrificed. To have continuity with respect to the data, we require that the position of the inserted knot tends to one of the end points of the interval as the data passes from requiring two cubic pieces to requiring only one. In addition, the slope d* at the inserted knot must tend to the slope at the corresponding end point.
21
To get a more visually pleasing curve, we allow A~ and A2 to move "away" from A, to avoid a flat curve. First we investigate the behaviour when h _< d2 -< ½ ( 3 A - d,). As d2 --~ ½(3A - d~), we note that the upper b o u n d for A~ (as given by Theorem 1 ) tends to A. At the same time, the range for A2 is bounded away from A and hence the knot x* tends to x2 (for any choice of &2 within the allowable range). Similarly d* --~ d2. This suggests that as d2 --~ ½(3A -- d~), we should choose A~ at its upper b o u n d so as to give continuity
8
7
6
5
4
3
3
2
= 1
0
i 2
i 3
J 4
? 5
•
-9 6
2
i 1
0
i 2
i 3
(a)
i 4
5
(b)
9,I 8
8
7
6
5
4
3
3
0
i 1
i 2
= 3
(c)
i 4
2
i
i
i
J
1
2
3
4
(d)
Fig. 7. (a) Osculatory Method; (b) Brodlie Formula; (c) Method A; and (d) Method B.
6
K. W. BRODLIE and S. BUTT
22
A similar a r g u m e n t applies when d2 > 3A - 2dr, then we are interested in continuity o f b e h a v i o u r as d2 -'-* 3A - 2 d l . The range for At is b o u n d e d away from A, while the lower b o u n d for A2 tends to A. Thus the k n o t tends to xt a n d continuity is assured. We choose
with respect to the data. However as d2 moves away from ½(3A - dr) a n d toward A, so we wish to choose At more toward the lower b o u n d o f the allowed range so as to avoid a fiat curve. T o further encourage this, A2 should be chosen at its u p p e r bound. T h u s we choose: A l = a U l + ( l -- a ) l t
A2 = u2
AI = II
d * = u3
A2 = flU2 + ( 1 -- fl)/2
d* = /3
where/3 = 1 - 3A - 2 d l / d 2 . Since there is no upper b o u n d for (/2 so we have constructed p a r a m e t e r fl ~ [0, 1] depending u p o n
where ut, u2, u3 are m a x i m u m o f the allowed ranges for At, A2, d* ranges respectively, a n d Ii is the minim u m of the range for At. The p a r a m e t e r a varies from 0 to 1 as d2 goes from A to ½(3A - d t ) a n d could be d e t e r m i n e d by c~ = 2(d2 - A ) / A - dr.
D = d2 - (3A - 2 d , )
20
20 f
17.5
17.5
15
15I
12.5
12.5
10
10
7.5
7.5
5
5
2.5
2.5
0¢
1
2
3
4
5
6
0
~
~ 1
2
(a)
,
t
3
4
5
* 6
(b)
20 17.5 15 12.5 10 7.5 5 2.5
2.5
1
2
3 (C)
4
5
5
0
1
i 2
i 3
J 4 (d)
Fig. 8. (a) Osculatory Method; (b) Brodlie Formula; (c) Method A; and (d) Method B.
i 5
J 6
Preserving convexity such that as D -+ O, [3 ~ 0 i.e., when the data approaches from requiring insertion of a knot to allowing a single convex cubic i.e., d2 -+ 3A -- 2d, from above, m i n i m u m values of Ai and A2 are attained, the inserted knot tends to x~ and d* converges to the slope at xt. This leads to the following method. Method B
In an interval where the cubic polynomial is nonconvex and d2 < ~ (3A - d t ), insert an extra knot with position of knot and the slope at knot determined by AI = aul + (I - a)l,
A2
=
u2
d*
=
bl 3
where a = 2(d2 - A ) / A -- d,. However if d2 > 3A - 2d~ then A, = l,
A2 = flu2 + ( l - f l ) h
d* =13
where fl = 1 - 3A - 2dL/d2. The resulting piecewise cubic polynomial is convex and the process is not sensitive to the change from a single cubic polynomial to a pair of cubics. Figure 6 shows the curve generated by Method B for our worked example, taking once again dE = 1.449 and 1.451. Two curves lie almost on top of one another. 6. NUMERICAL EXAMPLES
In this section we consider two example data sets and compare the results of the methods described in the previous sections with the osculatory method and the Brodlie formula (which preserves monotonicity but not convexity). For both Methods A and B, slopes at the nodes are estimated by means of the Brodlie formula. The first data set is given by
x0112+tnl+ r+ y
9
4
3
2.4
2.2
2.15
2.1
Results are shown in Fig. 7. The spurious inflection point in Fig. 7a introduced by the osculatory method is undesirable--the curve is known to be everywhere concave and free from inflection points. The interpolant (Fig. 7b) produced by the Brodlie formula is nonconvex as well. Both Methods A and B (Fig. 7c and 7d respectively) produce the desired curves. The second data set is given by xt01l y 0
0.5
2 13 14 I5 1.1 2.0 3.5 5.1
6 20
Results are shown in Fig. 8. Here the rather loose type of curve in Fig. 8a is generated by the osculatory method. Once again, interpolant Fig. 8b generated by Brodlie formula loses convexity. However, both Methods A and B (Fig. 8c and 8d respectively) produce
23
curves with the expected shapes by "repairing" the piece of the curve which loses convexity. 7. CONCLUSION & EXTENSIONS This paper has shown how a convex interpolant can be constructed, using only cubic polynomial pieces (results for polynomials of arbitrary degree will be presented in a subsequent paper). The slope at each data point is evaluated by some conventional method (such as the Brodlie formula). A simple test determines whether the cubic Hermite interpolant in an interval is convex. If not, a pair of cubic pieces can be used, and this paper has given conditions on the position of the knot and the slope of the interpolant at the knot. Two possible methods have been suggested: one which maintains second derivative continuity at the knot: the other which yields an interpolant that is continuous with respect to the data. The second method gives visually pleasing curves. Note that the method is also applicable to the case where derivatives are given explicitly rather than estimated. This is the case, for example, when plotting the solution components of ordinary differential equations[4 ]. REFERENCES
1. T. G. Ackland, On osculatory interpolation, J. Inst. Actuar. 49, 369-375. 2. J. H. Ahlberg, E. N. Nilson, and J. L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York and London, 12 (19XX). 3. J. Butland, Method for interpolating reasonable-shaped curves through any data, Proceeding of computer graphics, 80, Online Publications, 409-422 (1980). 4. R. W. Brankin and I. Gladwell, Shape preserving local interpolation for plotting solution of ODEs, 1MA Journal of Numerical Analysis. 9, 555-566 (1989). 5. K. W. Brodlie, A review of methods for curve and fimction drawing, K. W. Brodlie (Ed.), Academic Press, New York and London, 1-37 (19XX). 6. K. W. Brodlie, Methods for drawing curves, In Fundamental Algorithms for Computer Graphics, R. A. Earnshaw ( Ed. ), NATO ASI Series F17, Springer-Verlag, Berlin, 303-324 (1985). 7. R. L. Dougherty, R. L. Edelman, and J. M. Hyman, Positivity, monotonicity, convexity preserving cubic and quintic Hermite interpolation, Los Alamos National Lab-
oratory, LA-UR-85-2877 (19XX). 8. F. N. Fritsch and R. E. Carlson, Monotone piecewise cubic interpolation, S I A M Z Numer. Anal 17, 238-246 (19XX). 9. T. N. T. Goodman and K. Unsworth, Shape preserving interpolation by curvature continuous parametric curves, Computer Aided Geometric Design 5, 323-340 ( 1988 ). 10. J. A. Gregory and R. Delbourgo, Shape preserving interpolation, SIAM ,L Stat. Sci. Comp. 6, 967-976. 11. NAG Graphics Library Manual, Mark 3, NAG Ltd., Oxford, UK (1989). 12. E. Neuman, Uniform approximation by some Hermite interpolating splines, J. Comp. andAppl. Math. 4( 1), 79 (1978). 13. L. L. Schumaker, On shape preserving quadratic spline interpolation, SIAM J. Num. Anal. 20(4), ( 1983 ).
0097-8493/91 $3.00+ ,00 ~' 1991 Pergamon Press plc
Printed in Great Britain.
Technical Section
P R E S E R V I N G C O N V E X I T Y U S I N G PIECEWISE CUBIC I N T E R P O L A T I O N K. W. BRODLIE and S. BUTT School of Computer Studies, University of Leeds, Leeds, UK Abstract--This paper considers the problem of drawing a smooth cubic through a set of data points (x, i = 1, 2. . . . . N, where the y-values are dependent on the x-values--a common problem in scientific visualisation. A typical approach is to estimate the slope of the curve at each data point and construct a pieeewise cubic interpolant which is then easy to plot. An additional requirement is often to preserve the inherent shape of the data. While monotonicity can be preserved by suitable slope selection, it is known that the same cannot be achieved in general for convexity. However, this paper shows that by allowing (if necessary) two cubic pieces in some data intervals rather than one, convexity can always be preserved. Two methods are presented and illustrated with example data.
Yr),
I. I N T R O D U C T I O N
A fundamental problem in computer graphics is the drawing of a smooth curve through a set of data points (xr, Yr), i = 1, 2 . . . . . N. In many applications, particularly in scientific visualisation, the y-values are dependent on the x-values, and it is natural to seek an interpolant of the form:
y =y(x) which matches the given data. An attractive class of interpolating functions are piecewise cubic polynomials which are easy to construct and easy to plot. A typical procedure is to estimate the slope dr of the curve at each data point by some method; this then uniquely defines the cubic piece between each pair of data points. Many slope estimation methods have been suggested--a simple one is the osculatory, or cubic Bessel[l], method in which the slope at a data point is taken as the slope of a quadratic passing through the data point and its neighbour on either side. A full review of different methods is given in Brodlie[5, 6]. The methods all have first derivative continuity at the data points by construction. Attention has focussed recently on methods that preserve some shape inherent in the data. For example, it may be known that the function underlying the data is monotonic, and it is desirable to impose monotonicity as a constraint on the interpolant. The monotonicity case was explored by Fritsch and Carlson [ 8 ] who proved the following result: suppose the underlying function is monotonically increasing (i.e., Y l <- Y 2 " • • -< YN), then the cubic interpolant in each interval [xr, x~+~] is monotonically increasing if the slopes at the end-points satisfy
the above inequalities. This means that given any monotonic data it is possible to construct a monotonic piecewise cubic interpolant. The Brodlie formula is used in the N A G Graphical Library [ 11 ] curve drawing routines, and has proved visually acceptable. Another shape property of importance is convexity which is the subject of this paper. Here the situation is less satisfactory. In contrast to the monotonic case, there exist convex data sets (i.e., data from an underlying convex function) for which it is not possible to construct a convex piecewise cubic interpolant. This is discussed further in Section 2, together with various solutions to the problem generally involving a different class of interpolant to handle the impossible cases. In Section 3 we show that convexity can be achieved with piecewise cubics; the impossible cases are handled by using two cubic pieces in an interval between data points. In Sections 4 and 5 two specific methods are derived: one ensures second derivative continuity at the join between the two cubic pieces; the other ensures a smooth transition from the case where a single cubic suffices to the case where two pieces are required. Section 6 shows the methods applied to a number of test datasets that have appeared in the literature. 2. PRESERVING CONVEXITY IN INTERPOLATION Suppose the data points are realisations of some underlying convex function. Thus Ar+l >-Ai
i = 1,2 . . . . .
N-2.
(2.1)
Suppose we assign slopes dr to the data points. In the interval [xr, xr+l ], the cubic polynomial interpolant is convex [ 12 ] if and only if
di < 3Ai
2d; + d,+~ _< 3A; _< d; + 2di+~.
(2.2)
dr+l -< 3A,,
Handscomb (Brodlie[6]) points out that this implies, for five consecutive data points,
where Ai = (Y;+I - Y i ) l ( x r + l - xr). Butland[3] and Brodlie[6] give simple formulas for dr in terms of(xr, Yr) and its two neighbours that satisfy
Ai+ 2 --
15
Ai- I ~
1 . 5 ( A i + I -- A i )
(2.3)
K. W. BRODLIEand S. BUTT
16
which is a condition on the data--thus there will exist convex data for which there is no piecewise cubic interpolant. Therefore, there is no incentive to seek a "magic" formula that will generate slopes di always satisfying the convexity condition (2.2), as was possible in the monotonicity case. Various strategies for preserving convexity have been suggested in the literature. Schumaker [ 13 ] shows how quadratic spline interpolation can preserve convexity: slopes di are estimated by any method one chooses, with the restriction that d~+l > dj. Then in each interval [x~, xi+t ], a quadratic spline interpolant is calculated, having at most one internal knot. Schumaker gives conditions on the position of that knot that ensure that the slope d* at the knot satisfies:
di+j >- d* > di. The resulting quadratic spline is necessarily convex. The drawback of the method is the fact that it has lower order than cubic interpolation methods. Another approach is to use piecewise rational functions. Gregory and Delbourgo[10] and Brankin and Gladwell [ 4 ] both suggest the use of piecewise rational cubics. The drawback of this approach is the additional complexity of plotting a piecewise rational cubic. Indeed Brankin and Gladwell only recommend the rational cubic when the simple cubic fails to be convex. Still another approach is to use parametric curves: Goodman and Unsworth [9] describe an interpolation scheme using parametric cubic polynomials that can generate single-valued curves that preserve convexity. However, a disadvantage of using parametric curves for the single-valued case is that it is relatively difficult to evaluate y for a given value ofx. 3. PRESERVING CONVEXITY BY USING T W O CUBIC PIECES
xl = 0 ,
y~ = 0 ,
x2=l,
y2 = 1,
with estimated slope at xt as dj = 0.1. Then, from the convexity condition (2.2), we see that the cubic Hermite interpolant is convex provided 1.45 < d2 < 2.8. For values of d2 outside this range, the cubic Hermite is nonconvex, and Fig. 1 illustrates the situation for d2 =5.
Theorem 1 Let (x~, yj ), (x2, Y2) be two data points, and let d~, d2 be the (estimated) slopes at these points. Suppose the data is convex i.e., dt < a < d 2 but the associated cubic polynomial interpolant is nonconvex i.e., d2 ~ [ t ( 3 A - d l ) , 3A -- 2d,] where A = Y2 - Yl/X2 Xl. Then we can insert a knot at x* E [xl, x2] and associate with that knot a value y* and slope d*, such that the piecewise cubic interpolant fitted over the two intervals [xl, x* ], and [x*, x2] is convex. Moreover the range of knot positions and associated values and slopes is given by: -
-
Suppose the data are convex, and we use a slope estimation method which ensures 1
di+l>di
i = 1,2 . . . . . N - 1 . 0.8
(Note that every common method will do this.) If condition (2.2) is satisfied, then the cubic interpolant in [x~, x~+l] will preserve convexity. So we need only be concerned with the case where (2.2) is not satisfied. The resulting cubic in this case is nonconvex, and we need some way of incorporating an extra degree of flexibility. We would prefer to do this while still using piecewise cubics. Dougherty, Edelman, and Hyman[7] suggest the relaxing of first derivative continuity at the data points, but this could yield a curve that is not visually acceptable. Our suggestion is to create the flexibility by the insertion of an extra knot between [xi, x~+~]. At this stage we simplify the notation by considering only the interval [x~, x2]. It will be useful also to have a specific example to consider in parallel with the theoretical development. We consider the case:
0.6
0.4
0.2
-012
O4
06 /
O8
-0.2
Fig. 1. Nonconvex cubic hermite interpolant.
t
Preserving convexity A, ~ [ m a x { d ~ , A -
](d2-
17
and so we extend ( 3.1 ) to
d,)},
m i n { A , d 2 - ~ ( d 2 - d~)}]
d~ < A, < min { A, d2 - ~ ( d2 - d, ) } .
(3.11)
A2 E [ m a x { A , A, + ~ ( d 2 - d , ) } , min {d2, A~ + ](d2 - d , ) } l d* E [max{ ½(3A~ - d,), 3 A 2 - 2d2},
( N o t e d~ < d2 - ~(d2 - d j ) , so (3.11) is c o n s i s t e n t . )
Next consider inequality (3.9). This can be written as:
m i n { 3 A , - 2d~, ½(3A2 - d2)}] where A, = (y* - y l ) / ( x * - x l ) , A2 = (Y2 - y * ) / ( x 2 -- X*).
A2 < A, + ](d2 - d~) and so we extend (3.10) to
Proof First, A, and A2 must satisfy: m a x { A , A, + ~ ( d 2 - d,)} < A2 d, < A, < A
(3.1)
A < A2 < d2
(3.2)
< min{d2, A, + ] ( d 2 -
dr)}.
(3.12)
This requires for consistency that Next d* must satisfy: A-
(3.3)
½(3AI -- d l ) < d* < 3A, - 2d,
(3.4)
½(3A2 - d*) -< d2 -< 3A2 - 2d*.
(3.5)
or
Al>A-~(d2-dl)
and so we extend ( 3.11 ) to max{ dl, A - ~(d2 - d l ) } , -< AI <-min{A, dz-~(d2-d,)} If inequalities ( 3 . 1 ) - ( 3.5 ) are satisfied, then the piecewise cubic interpolant is convex. Inequality (3.5) can be written as: 3A2 - 2d2 -< d* < ½(3A2 - d2)
(3.6)
and inequalities (3.3), (3.4), and (3.6) can be combined into the single inequality
which is easily shown to be consistent. Hence there exists a A~ satisfying (3.13), a A2 satisfying ( 3.12 ) and d* satisfying (3.7), and the theorem is proved. Note that the position of the knot (x*, y*) will be given by
(3.7)
using the fact that AI < ½(3Al -- d,) and A2 > ½(3Az - d2) which is true ifAl and A2 satisfy (3.1) and (3.2) respectively. Again if (3.1) and (3.2) are true, there will be a value of d* satisfying (3.7) provided that: ½(3A, - d,) < ½(3A2 - dE)
(3.8)
3A2 - 2d2 < 3Al - 2dl.
(3.9)
1
x*
- {(x2A2 - x l A l ) - - (Y2 -- Yl)} A2 - AI
y*
- { AiA2(x2 - Xl) - (A, y2 - A2yl)}. A2 - A3 (3.15)
max{ ½(3A, - d,), 3A2 - 2d2} ~ d* -< min{ 3A, - 2d,, ½(3A2 - d2)}
(3.13)
(3.14)
1
We can illustrate this theorem in terms of our worked example. Consider a slope d2 = 5.0 which causes the convexity test to fail. The theorem states that a knot can be inserted with associated slope d*, such that
+ )min(10 ,O;)]
Consider inequality (3.8) first. This can be written as i.e.,
A~ E5 [0.1, 1.01,
A2 ~ A, + ](d2 - d,) A2 ~ [ m a x ( 1.0, A, + - ~ ) , min(5.0, A, + ? ) ]
and we can extend (3.2) to m a x { A , A, + ](dz - d,)} -< Az ~ d2.
(3.10)
This requires for consistency that A,+ ](d2-d,)
or
A, < d 2 - ] ( d 2 - d,)
i.e.,
A2E[A,+-~,A,+?],
d* E [max{ ½(3At - 0.1), 3A2 - 10.0}, min{3A~ - 0.2, ½(3A2 - 5.0)}].
K. W. BRODLIEand S. BUTt
18
In Fig. 2, we choose A, = 0.5, A2 = 2.5, and d* = 1.0 and draw the resulting piecewise cubic which is convex. A knot is inserted as shown at (0.75, 0.375). Theorem 1 shows there is a range of possible knots and associated slopes. The following corollary shows that a simple choice of A,, A2 always lies within the range given in the theorem. Corollary l: Let A, = a A + (1 - a ) d i ,
A2 = a A + (1 - a)d2.
we require aA+(l
-a)d,
> A-
2(d2-dt)
i.e.,
2(d2-dl)>(1-a)(A-d,)
which is true for a E [ }, 2 ] since d2 > A. To show A , <--d 2 - -
]1( d 2 - d , )
we require
If a E [~, 2 ], then A1, A2 lie within the range given in Theorem 1 for all dr, d2, A, and the corresponding range for d* becomes:
aA + (1 -- a ) d , < d2 - ](d2 - d,) i.e.,
](d2-d,)>a(A-d,)
which is true for a E [0, ] ] . Thus
d* ~ [ m a x { d , + 3a(A - dj ), d2 - 3a( d2 - A)},
A, e [max{d,, A - ](d2- d,)},
min {d, + 3a(A - d,), d2 - 3a(d2 - A)}].
min{ A, d2 - ](d2 -
d,)}l
Proof." We have to show that
f o r a E [I, 21. Next we have to show that A, E [ m a x { d , , A -
2(d2 - d,)}, A2 E [max{A, A, + l ( d 2 -
min{ A, d2 - l ( d 2 - d , ) } ] .
d,)}, min {d2, Al](d2 - dj)}]
The results:
The results A, >_ d~
and
A, -~ A
A2> A
are trivial. To show A,>A-
and
A2
are trivial. To show A2 > A, + l ( d 2 - d , )
2(d2-d,) we require aA+(1
-- o~)d2 > aA + (1 - a ) d , i.e..
0.8
+ ~(d2-d~)
(1 - a ) ( d 2 - d , ) ~
~(d2-dt)
which is true for a E [ ~, 1]. To show A 2 <_ A, + 2(d2 - d , )
0.6
we require 0.4
otA + (1 - a)d2 -< aA + (1 -- a ) d , + ~(d2 - dr) i.e.,
(1-a)(d2-d,)_<2(d2-d,)
0.2
which is true for a E [0, 2]. Thus °
-0.2
0.2
0.4
0.6
0.8
1
A2 ~ [max{A, A, + -~(d2-
d,)}, min{d2, A,~(d2 - d,)}]
-0.2
Fig. 2. Convexity preserved by inserting a knot.
for a E [I, 21.
Preserving convexity The range for d* follows by simple substitution. Note: The extreme choices a = ] and a = 2 lead to particularly simple ranges for d * , viz..
a = ~ --,- d* ~ [ m a x { ' i ( d l + A), A}, min{A, ½(A+d2)}] a = ]--~d*@[max{A,
i.e.,
d* = A
2A-d2},
19
C o r o l l a r y 2: Let Al = a A + ( l - a ) d l ,
A2 = a A + ( 1
- a)d2
as in Corollary 1, with cr E [ -~, -] ]. Then the required value of d* for second derivative continuity is d* = [(d, + d2) - A] + 3a[A -- ½(d, + d2)]
min{2A--d1,A}]
i.e.,
d * = A.
4. A M E T H O D ENSURING SECOND DERIVATIVE CONTINUITY AT KNOT
To obtain a particular algorithm that will maintain convexity, we need to make a specific choice of A~, A2, and d* satisfying Theorem 1. It is convenient to work with the simpler formulas of Corollary 1. An attractive property is that the two cubic pieces be joined smoothly at the knot x * i.e., we seek second derivative continuity at x*, in which case the curve is a clamped spline. We have the following theorem.
We now explore for which values of a the value of d* given in Corollary 2 lies with the valid range of d* as given by Corollary 1. We have the following results. Theorem
3
Let Aj, A2 be as in Corollary 1. Then the value of d* giving second derivative continuity as identified in Corollary 2 lies within the valid range for all possible dl, d2 and A if and only if a = -~. P r o t f : Suppose a = ~. Then d* = A and it has already been noted after Corollary 1 that this is valid. Suppose a = ~ - e, for some { > 0. Then
Theorem 2
Consider a piecewise cubic c ( x ) constructed over the interval [x~, x2] so that
d* = a -
3~{a-
½(d, + d2)}.
Corollary 1 requires, among other things, c ( x l ) = )'j,
c(x2) = Y2,
c'(xl)
= d~,
c ' ( x 2 ) = d2
and with a knot at x* where x~ < x* < x2. Then the piecewise cubic has a continuous second derivative at x* if and only if
d*>_d2-3a(d2-A)=A-(1-30(d2-A). This requires: 3~(d2 - d l ) < d2 - A.
d ~
--
2 (x: - x, )
{3(hA 2 + k A , ) - ( h d 2 + k d j ) }
where h = x* - xl, k = x2 - x*. P r o o f Let c ~ ( x ) be the cubic piece over [x,, x * ] and c2(x) the cubic piece over I x * , x2]. Then it is easy to show (see Ahlberg, Nilson, and Walsh[2]) that 2dj 4d* c'~{x* ) = T + T 2(x
) -
4d*
60'2- y*)
k
k
k2
Method A
In any interval where the cubic polynomial is nonconvex, insert an extra knot with position determined by
6(y* - Yl) h~
2d2
Given any ~ > 0, we can always choose d~ < A < d2 such that this inequality fails. Hence the theorem is proved. This leads to the following method.
A, = 2A + ~d,
(4.1)
A2 = -~A + ~d2
(4.2)
Thus for c'f(x*) = c~(x*), we require
t~- + 4d*
+
-~ k _ 6(y* - yl) 4 6(y2 - y * ) h2
i.e., d * - - -
1
2(h + k)
and assign a slope d* = A at the knot. The resulting piecewise cubic polynomial is convex and is a cubic spline. Substituting (4.1) and (4.2) into (3.14) and (3.15) give the knot position as
k2
1
{ 3(hA2 + kay) - ( h d 2 + k d , ) } .
Theorem 2 tells us the following: once we have selected a knot x* and associated y*, we must choose a specific value of d* to obtain second derivative continuity.
X*
- {(x2d2 - x , d , ) d2 - dl
- (y2 - yl)}
y* = - [{d, d2 + 2A(d, + d2 - a ) } ( x 2 - x,) d2 - dl - 3 ( d , y2 - d2y,)].
20
K.W.
BRODLIEandS.
1
BUTT
1
0.8
0.8
0.6
0.6
0.4
0.4 0.2
0.2 i
0.2
-0.2
0.4
0.6
0.8
1
0.2
-0.i
Fig. 3. Method A--preserving 2nd order continuity at knot.
Figure 3 shows the effect of applying this method to our worked example. A knot is inserted as shown at (x*, y * ) = (0.82, 0.57) with slope d* = 1.0. 5. A METHOD INSENSITIVE TO CHANGES IN DATA Method A suffers from two difficulties. Firstly, it is a desirable property of any curve drawing method that the generated curve should not be sensitive to small changes in the data. Unfortunately this is not the case with Method A as the data moves from allowing a single convex cubic to requiring the addition of an internal knot. This is illustrated in Fig. 4, where the
0.4
0.6
0.8
1
Fig. 5. Always inserting a knot.
curves for the cases dE = 1.451 and d2 = 1.449 are shown alongside each other. This difficulty can be overcome by always inserting a knot, even when a single cubic will suffice. That is, we calculate Al and A2 by (4.1) and (4.2), and take d* = A, and use a pair of cubic pieces whatever the data. Since the expressions for At, A2 and d* are continuous with respect to the data, the resulting interpolant will likewise be continuous with respect to the data. This is illustrated in Fig. 5, where the curves for the cases d2 = 1.451 and d2 = 1.449 are seen to lie almost on top of one another. Contrast this with Fig. 4. However, there is a second difficulty, and this is the visual quality of the curve. The choice a = ] means that Aj and A2 are both relatively close to A, and so 1
0.8
0.8
0.6
/
0.4
,,," 06
0.4
0.2
o.~
o.',
o18
o'.8
;
0.2
C l a m p e d cubic spline Single cubic Fig. 4. Sensitivity to change of data.
OI
0.2
0.4
0.6
0.8
Fig. 6. Insensitive to change of data.
1
Preserving convexity the inserted knot is close to the chord joining the two data points and a "flat" curve results. In this section, therefore, we discuss an alternative method which is not sensitive to small changes in the data, and which aims to give a more visually pleasing curve. Second derivative continuity at the knot is sacrificed. To have continuity with respect to the data, we require that the position of the inserted knot tends to one of the end points of the interval as the data passes from requiring two cubic pieces to requiring only one. In addition, the slope d* at the inserted knot must tend to the slope at the corresponding end point.
21
To get a more visually pleasing curve, we allow A~ and A2 to move "away" from A, to avoid a flat curve. First we investigate the behaviour when h _< d2 -< ½ ( 3 A - d,). As d2 --~ ½(3A - d~), we note that the upper b o u n d for A~ (as given by Theorem 1 ) tends to A. At the same time, the range for A2 is bounded away from A and hence the knot x* tends to x2 (for any choice of &2 within the allowable range). Similarly d* --~ d2. This suggests that as d2 --~ ½(3A -- d~), we should choose A~ at its upper b o u n d so as to give continuity
8
7
6
5
4
3
3
2
= 1
0
i 2
i 3
J 4
? 5
•
-9 6
2
i 1
0
i 2
i 3
(a)
i 4
5
(b)
9,I 8
8
7
6
5
4
3
3
0
i 1
i 2
= 3
(c)
i 4
2
i
i
i
J
1
2
3
4
(d)
Fig. 7. (a) Osculatory Method; (b) Brodlie Formula; (c) Method A; and (d) Method B.
6
K. W. BRODLIE and S. BUTT
22
A similar a r g u m e n t applies when d2 > 3A - 2dr, then we are interested in continuity o f b e h a v i o u r as d2 -'-* 3A - 2 d l . The range for At is b o u n d e d away from A, while the lower b o u n d for A2 tends to A. Thus the k n o t tends to xt a n d continuity is assured. We choose
with respect to the data. However as d2 moves away from ½(3A - dr) a n d toward A, so we wish to choose At more toward the lower b o u n d o f the allowed range so as to avoid a fiat curve. T o further encourage this, A2 should be chosen at its u p p e r bound. T h u s we choose: A l = a U l + ( l -- a ) l t
A2 = u2
AI = II
d * = u3
A2 = flU2 + ( 1 -- fl)/2
d* = /3
where/3 = 1 - 3A - 2 d l / d 2 . Since there is no upper b o u n d for (/2 so we have constructed p a r a m e t e r fl ~ [0, 1] depending u p o n
where ut, u2, u3 are m a x i m u m o f the allowed ranges for At, A2, d* ranges respectively, a n d Ii is the minim u m of the range for At. The p a r a m e t e r a varies from 0 to 1 as d2 goes from A to ½(3A - d t ) a n d could be d e t e r m i n e d by c~ = 2(d2 - A ) / A - dr.
D = d2 - (3A - 2 d , )
20
20 f
17.5
17.5
15
15I
12.5
12.5
10
10
7.5
7.5
5
5
2.5
2.5
0¢
1
2
3
4
5
6
0
~
~ 1
2
(a)
,
t
3
4
5
* 6
(b)
20 17.5 15 12.5 10 7.5 5 2.5
2.5
1
2
3 (C)
4
5
5
0
1
i 2
i 3
J 4 (d)
Fig. 8. (a) Osculatory Method; (b) Brodlie Formula; (c) Method A; and (d) Method B.
i 5
J 6
Preserving convexity such that as D -+ O, [3 ~ 0 i.e., when the data approaches from requiring insertion of a knot to allowing a single convex cubic i.e., d2 -+ 3A -- 2d, from above, m i n i m u m values of Ai and A2 are attained, the inserted knot tends to x~ and d* converges to the slope at xt. This leads to the following method. Method B
In an interval where the cubic polynomial is nonconvex and d2 < ~ (3A - d t ), insert an extra knot with position of knot and the slope at knot determined by AI = aul + (I - a)l,
A2
=
u2
d*
=
bl 3
where a = 2(d2 - A ) / A -- d,. However if d2 > 3A - 2d~ then A, = l,
A2 = flu2 + ( l - f l ) h
d* =13
where fl = 1 - 3A - 2dL/d2. The resulting piecewise cubic polynomial is convex and the process is not sensitive to the change from a single cubic polynomial to a pair of cubics. Figure 6 shows the curve generated by Method B for our worked example, taking once again dE = 1.449 and 1.451. Two curves lie almost on top of one another. 6. NUMERICAL EXAMPLES
In this section we consider two example data sets and compare the results of the methods described in the previous sections with the osculatory method and the Brodlie formula (which preserves monotonicity but not convexity). For both Methods A and B, slopes at the nodes are estimated by means of the Brodlie formula. The first data set is given by
x0112+tnl+ r+ y
9
4
3
2.4
2.2
2.15
2.1
Results are shown in Fig. 7. The spurious inflection point in Fig. 7a introduced by the osculatory method is undesirable--the curve is known to be everywhere concave and free from inflection points. The interpolant (Fig. 7b) produced by the Brodlie formula is nonconvex as well. Both Methods A and B (Fig. 7c and 7d respectively) produce the desired curves. The second data set is given by xt01l y 0
0.5
2 13 14 I5 1.1 2.0 3.5 5.1
6 20
Results are shown in Fig. 8. Here the rather loose type of curve in Fig. 8a is generated by the osculatory method. Once again, interpolant Fig. 8b generated by Brodlie formula loses convexity. However, both Methods A and B (Fig. 8c and 8d respectively) produce
23
curves with the expected shapes by "repairing" the piece of the curve which loses convexity. 7. CONCLUSION & EXTENSIONS This paper has shown how a convex interpolant can be constructed, using only cubic polynomial pieces (results for polynomials of arbitrary degree will be presented in a subsequent paper). The slope at each data point is evaluated by some conventional method (such as the Brodlie formula). A simple test determines whether the cubic Hermite interpolant in an interval is convex. If not, a pair of cubic pieces can be used, and this paper has given conditions on the position of the knot and the slope of the interpolant at the knot. Two possible methods have been suggested: one which maintains second derivative continuity at the knot: the other which yields an interpolant that is continuous with respect to the data. The second method gives visually pleasing curves. Note that the method is also applicable to the case where derivatives are given explicitly rather than estimated. This is the case, for example, when plotting the solution components of ordinary differential equations[4 ]. REFERENCES
1. T. G. Ackland, On osculatory interpolation, J. Inst. Actuar. 49, 369-375. 2. J. H. Ahlberg, E. N. Nilson, and J. L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York and London, 12 (19XX). 3. J. Butland, Method for interpolating reasonable-shaped curves through any data, Proceeding of computer graphics, 80, Online Publications, 409-422 (1980). 4. R. W. Brankin and I. Gladwell, Shape preserving local interpolation for plotting solution of ODEs, 1MA Journal of Numerical Analysis. 9, 555-566 (1989). 5. K. W. Brodlie, A review of methods for curve and fimction drawing, K. W. Brodlie (Ed.), Academic Press, New York and London, 1-37 (19XX). 6. K. W. Brodlie, Methods for drawing curves, In Fundamental Algorithms for Computer Graphics, R. A. Earnshaw ( Ed. ), NATO ASI Series F17, Springer-Verlag, Berlin, 303-324 (1985). 7. R. L. Dougherty, R. L. Edelman, and J. M. Hyman, Positivity, monotonicity, convexity preserving cubic and quintic Hermite interpolation, Los Alamos National Lab-
oratory, LA-UR-85-2877 (19XX). 8. F. N. Fritsch and R. E. Carlson, Monotone piecewise cubic interpolation, S I A M Z Numer. Anal 17, 238-246 (19XX). 9. T. N. T. Goodman and K. Unsworth, Shape preserving interpolation by curvature continuous parametric curves, Computer Aided Geometric Design 5, 323-340 ( 1988 ). 10. J. A. Gregory and R. Delbourgo, Shape preserving interpolation, SIAM ,L Stat. Sci. Comp. 6, 967-976. 11. NAG Graphics Library Manual, Mark 3, NAG Ltd., Oxford, UK (1989). 12. E. Neuman, Uniform approximation by some Hermite interpolating splines, J. Comp. andAppl. Math. 4( 1), 79 (1978). 13. L. L. Schumaker, On shape preserving quadratic spline interpolation, SIAM J. Num. Anal. 20(4), ( 1983 ).