Monotone linear rational spline interpolation

Computer

Aided

Geometric

313

Design 9 (1992) 313-319

North-Holland

COMAID

289

Monotone linear rational spline interpolation Richard D. Fuhr and Michael Kallay Electronic Data Systems, 13555 SE 36th St. Suite 300, Bellcue, WA 98006, USA Received

February

Revised April

1991

1992

Abstract Fuhr, R.D. 9 (1992)

and M. Kallay, Monotone

linear rational

spline interpolation,

Computer

Aided

Geometric

Design

with a C’

monotone

rational

313-319.

A method

for interpolating

B-spline of degree ing B-spline

a monotone

1 is presented,

curves in Geometric

Keywords. Interpolation;

data sequence

(values and derivatives)

along with some test results. This method

is tailor-made

for reparameteriz-

Modeling.

monotonic;

reparameterization;

splines

1. Introduction In geometric modeling, it is often useful to be able to reparameterize curves; e.g., when constructing a surface from these curves. Let C(t) be a parametric curve, defined for a < t < b. We can construct a reparameterized version c of C by finding a monotone function f(s) with f: [(u, PI + [a, bl such that d(s) = C(f(s)). The application will determine what smoothness requirements should be imposed on f. One application of reparameterization is to make a G’ curve C’. For a simple example, suppose that C: [0, 11+ R3 and suppose that C’(t) is continuous for all t E 10, l] except for one a E [0, l] where C’(a - ) = kC’(a + 1 for some k f 0. Our task is to find a monotone function f(s) with f: [O, 11--f [O, 11 such that f’) is a reparameterized version of C that is CL everywhere. Another use for reparameterization is in the construction of ruled surfaces where ‘matched points’ are specified. Given rail curves C,(S) and C,(t), we identify a list of pairs of parameters ((si, t,): i = 0,. . . , n) that are used to specify pairs of matched points {(C,(s,), C&t,))}; i.e., for each i, a ruling line in the surface connects the points C,(si) and C,(ti). Here, we assume that o=s,
..I


and

O=t,
e.1
We also assume that C, and C, are parameterized on [0, 11. To reparameterize C, and C,, let ui = (si + ti)/2 for i = 0,. . . , n. We need to construct C’ monotone functions f,: [O, 11--) [O, 11 and fi: [O, 11--* [O, 11such that f,(ui> = si and fi = ti Correspondence

IO: M. Kallay, Electronic

0167-8396/92/$05.00

0

1992 - Elsevier

Data

Systems, 13555 SE 36th St. Suite 300, Bellvue, WA

Science Publishers

B.V. All rights reserved

98006,

USA.

314

R. D. Fuhr and M. Kallay / .Cionotone linear rational spline interpolation

for i = 0,. ..,n. We can use the reparameterized C,(fZ(u>> to build the ruled surface

curves d,(u) = C,(f,(u))

and Cz(~O =

F(u,c)=(l-L’)6,(rt)+cd,(u). We use the term rational B-spline to mean a vector-valued the form

(or scalar-valued) function F of

F( ~4)= E w,Pibi( U)/E Wlbi(U) where the 6, are a set of B-spline basis functions determined by a specified degree d and knot sequence ui. The Pi are vectors (or scalars) and the wi are positive scalars. For more details on rational B-splines see [Boehm et al. ‘841. We use the term single-span rational B-spline to mean a rational B-spline for which the underlying basis functions are the Bernstein polynomials. To construct a set of such basis functions of degree d on the domain [a, 61 the knot sequence would be a ,..., u,b ,..., b. -d+l

d+l

If the original curve C and the reparameterization function f are rational B-splines, then the reparameterized curve d = C 3 f is also a rational B-spline. The degree of d is the product of the degrees of C and f. To minimize the degree of 6, we seek to make the degree of f as low as possible; and we can’t do better than degree 1. To obtain the additional degrees of freedom necessary for C’ interpolation, we look at linear rutionuf B-splines (LRBS); i.e., rational B-splines of degree 1.

2. Previous work Monotone interpolating schemes for data with derivatives have been presented in [Fritsch & Carlson ‘SO],[Fritsch & Butland ‘841 and [Gregory & Delbourgo ‘821. The Fritsch-Carlson and Fritsch-Butland splines are cubic. Since the cubic Hermite polynomial is not necessarily monotone, even with increasing data and positive derivatives, the interpolating schemes in [Fritsch & Butland ‘841 and [Fritsch & Carlson ‘801 will modify the input derivatives if they violate the necessary condition for monotonicity. The newer Gregory-Delbourgo spline in [Gregory & Delbourgo ‘821 is quadratic rational. It will accommodate any set of derivatives whose signs are compatible with the monotonicity of the data. Our LRBS, like the Gregory-Delbourgo spline, will accept any set of monotonic data with compatible derivatives, but it cannot handle zero slope. This is no handicap for an interpolator used to reparameterize curves, because a vanishing derivative is undesirable in curve parameterization anyway. Compared to these schemes, our LRBS has the advantage that it doesn’t raise the degree of a curve when used for reparameterization. With the same data, the LRBS will require less storage space than the cubic and the rational quadratic monotone splines. When derivatives are not prescribed with the data, Gregory and Delbourgo have shown in [Gregory & Delbourgo ‘831 how to assign derivatives to the data points to obtain a C* interpolating quadratic rational spline. There is no way to do it with LRBS, because the second derivative of the LRBS cannot vanish, so it will not handle inflection points. An explicit expression for reparameterization with LRBS is given in [Lee& Lucian ‘911. Our focus here is not on curve reparameterization but rather on the C’ monotone interpolation with LRBS, which is required for effective C’ reparameterization.

R.D. Fuhr and M Kailay / Monotone linear rational splinr interpolation

315

3. The linear rational function We begin with a brief study of the behavior of the linear rational function (LRF) y(x)

= g+

=

W,Y,(b

-xl + WbYb(X - 0)

w,(b-x)

+Wb(X-u)

.

Here, we are specifying that f(x) is the numerator and g(x) is the denominator of the right-hand expression. This is the single span rational B-spline form with distinct knots u, b, coefficients y,, y, and weights wQ, wb (presumed positive). Note that f’ = wbyb - w,y,, and g’ = wb - w,. The first derivative of y is: y,=f’g-fg’

f’-yg’ Z----E g

g2 =

wb(yb -y)

Wbyb-W,y,-y(Wb-W,) g

+ w,(y -ya) LT

y is a convex combination

of y, and y,, so it lies between them. The sign of the numerator is therefore equal to the sign of (yb - y,), and y is strictly monotone throughout [a, b]. Since f(a) = w,y,, f(b) = why,, da> = wo, and g(b) = wb, the derivatives at the endpoints are: y’(u)

wb yb -ya = kb-a,

y’(b)

and

wa yb-y, = -~ wb b-u

’

Thus y’ never vanishes in the interval [a, 61, unless y is constant. Moreover, since f” = g” = 0, we have y” = -2g’y’/g, hence y” never vanishes on [a, 61, unless g’ = 0, which happens only when y is linear, because g’ = wb - w,. The linear rational B-spline (LRBS) y(x) with distinct knots x0,. . . , x,, coefficients Y,, . . . ,y,, and weights wo, . . . , w,, is strictly increasing if y, y, > . . . > y,. In any case, its first and second derivatives never vanish. Its one-sided first derivatives at the knot xi are: wi-I

Yi -Yi-1

wi

xi-xi_,

Y'(Xi -)

= -

y’(x;+)

wi+I = wi

,

(1)

-xi.

(2)

Yi+l -Yi xi+1

4. Interpolation with linear rational B-splines In general there is no linear rational function interpolating the values and derivatives ((a, y,, yi), (b, y,, y;)}, because we have three degrees of freedom against four constraints. But if we turn to a linear rational B-spline of two spans, then the number of free parameters will match the number of constraints. It turns out that these constraints can always be satisfied if the data is strictly monotone, i.e., positive derivatives and yb > y,,, or negative derivatives and y, < y,. For any h between 0 and 1 let c = (1 - h)u + Ab. The linear rational B-spline over the knot sequence (a, a, c, b, 6) is defined by three coefficients and three weights. The first and last coefficients must be y, and y, respectively, and we are free to choose any positive value for the first weight w,. Three conditions on the derivatives at a, b, and c will determine the

R.D. Fuhr and At. Kallay / Monotone linear rarronal spline inrerpolation

316

coefficient are:

y, and the other two weights wC and wb. The derivative interpolation

constraints

WCYc -Y,

yl=wa c-a

(3)

we yb-yc y;,=w, b-c.

(4)

and

The derivative continuity condition at c is

waYc-Y, --=--

wb yb

WC c-a

-yc (5)

wc b-c’

The system of equations (3)-W is nonlinear in the unknowns yC, w,, wb, but it can be solved analytically. Equation (5) can be solved for y,, after substituting b - c = Mb - a) and c-a=(l-AXb-a): Y,

=

(l-A)w,Y,+hWbYb (1 -A)w,+Aw,

(6)

’

Dividing (3) by (4) yields Y:,

wb

(yc

YL

w,

(yb-Yc)(C-a)’

-=-

-%T)(~

-c)

On the other hand, we have from (5)

hence

with the positive solution Wb

(7)

=

Finally write (3) and (4) as: wC(yC- y,) = w,yi(l - AXb - a) and wc(yb respectively, and add them up to obtain: w,=(Aw,y;+(l-A)w,yl,)-.

-Y,)

= wbYkA(b

- a>,

b-a yb-y,

The algorithm for constructing the spline is straightforward. wa is at our disposal, and wb is defined by (7), assuming monotone data yi/yL > 0. y, is defined by (61, again relying on the strict monotonicity of the data; its value will be between y, and yb. Finally (8) defines w,; it will be positive if the data is monotone. Monotonicity is thus necessary as well as sufficient: if the data is not monotone then the equations may have no solution at all, but if it is monotone then so will be the solution spline.

R.D. Fuhr and M. Kallay / Monotone linear rational spline interpolation

317

5. The choice of derivatives

If we are using an interpolating scheme that requires function values and derivatives, but we only have the data for function values, we need to specify suitable derivative information. Rules that associate derivatives with data points typically involve estimates based on slopes or function differences nearby. The reader may find a survey of various derivative-assignment schemes in Section 13 of [Boehm et al. ‘841. In all the schemes surveyed, the derivative assigned to xi is invariably a linear combination of nearby chord-slopes or data differences; the various schemes are evaluated based on their invariance under various geometric transformations-Euclidean, affine or projective. This approach makes sense in the context of geometric curves interpolating geometric points. However, a different consideration should determine the derivatives in the context of curve reparameterization. If one curve is reparameterized to suit another, for example. then it is desirable to make the process invariant under switching the curves. This invariance is not always possible; the inverse of a quadratic or cubic spline is not even rational. But any derivative-assignment scheme should aim to get close to this invariance: the derivative assigned to the (y,, xi) should be the reciprocal of that assigned to (xi, yi> and the derivative assigned to (xi, zi> should be the product of the derivatives assigned to (x1, yi> and (y,, zi>. This calls for assigning some mrlltiplicatice combination of chord-slopes (e.g., some weighted geometric mean), rather than a linear combination. In our scheme we chose y( to be the derivative of the unique LRF through the three data

0

10002

10.003 10005

LO.006 IO.008

Fig. 1. Akima’s

10.5

data.

I.5

50

60

85

318

points (xl_,,

R.D. Fuhr and M. Kallay / ,Wonotone linear rational spline interpolation

Y,_,

>, (xi, y,), (xi+ ,, y,+ I) for interior 0 < i < n. The formula turns out to be

very simple: y[ = ri- tSi/fi,

where ~,=(yi+,-Y~)/(~i+,-~i), and

ti=(Y,+1_Yi-,)/(xi+,-xi-,).

At the endpoints we chose yh to be the derivative of the LRF through (x,, ya), (x,, y,>, (x,, yz) and yi as that of the LRF through (x,_*, Y~_~), (x,-r, y,,_t), (x,,, Y,). This prescribes y; =~at,/s, and y; =~,_t~,,_,/s,_~.

6. An implementation We have implemented a simple C’ interpolating scheme based on this theory. Given the data {(xi, yi, yl): i = 0,. . . , n], we insert an additional knot at each span’s midpoint (A = +I. We set wa = 1, and then we use equations (6148) to compute the weights and the missing coefficients of the interpolating spline. The input may prescribe two distinct one-sided derivatives at any interior xi, to model a derivative discontinuity. If no derivatives (yi} are prescribed then we assign derivatives according to the scheme described above. To compare the quality of the interpolation, we tested our scheme on the data sets given in [Fritsch & Carlson ‘SO]. We had to cheat on the first 6 data points of Akima’s data [Akima ‘701, because they present zero slope, which no LRBS can handle. We tweaked those points to introduce a slope of 0.001. The result, which seems to be quite acceptable, is shown in Fig. 1. No cheating was necessary for the LLL (Lawrence Livermore Laboratories) radiochemical data [Fritsch & Carlson ‘801. Our interpolant is shown in Fig. 2.

1.99 0

a 09

8.19

2.76429 E-5

4.37498 E-2

a7 0.169La3

9.2 0.469428

LO 0.943740

I2

I5

20

0.998636

0.999919

a.999994

Fig. 2. LLL radiochemical data.

R.D. Fuhr and !M. Kal1a.v / ‘Lionorone linear rational splinr interpolarion

319

References Akima, H. (1970). A new method of interpolation and smooth curve fitting based on local procedures, J. Assoc. Comput. Mach. 17, 589-602. Boehm, W., Farin. G. and Kahmann, J. (19&t). A survey of curve and surface methods in CAGD. Computer Aided Geometric Design 1, l-60. Fritsch, F.N. and Butland, J. (19841, A method for constructing local monotone piecewise cubic interpolants. SIAM J. Sci. Stat. Comp. 5, 300-304. Fritsch, F.N. and Carlson, R.E. (19801, Monotone piecewise cubic interpolation, SIAM J. Numer. Anal. 17, 238-246. Gregory, J.A. and Delbourgo. R. (19821, Piecewise rational quadratic interpolation to monotonic data, IMA J. Numer. Anal. 2, 123-130. Gregory, J.A. and Delbourgo. R. (19831, Cz rational quadratic spline interpolation for monotonic data. IMA J. Numer. Anal. 3, 141-152. Lee, E.T.Y. and Lucian, M.L. (19911, Mobius reparameterization of rational B-splines, Computer Aided Geometric Design 8, 213-215.