Shape preserving approximation

Shape preserving approximation

435 Short Note Shape Preserving Approximation Jernej K o z a k Dept. of Mathematics and Mechanics, Edvard Kardelj Universtty, Jadranska 19, 61111 L/...

310KB Sizes 0 Downloads 54 Views

435

Short Note

Shape Preserving Approximation Jernej K o z a k Dept. of Mathematics and Mechanics, Edvard Kardelj Universtty, Jadranska 19, 61111 L/ubljana, Yugoslauia An algorithm that produces a shape preserving interpolatory hyperbolic spline~,is discussed here. It is based upon Schweikert's spline in tension, but with tension parameters determined in advance. These are chosen in such a way that the approximating function differs as little as possible from a cubic polynomial interpolant. By this approach the approximation power of the cubic spline is preserved wherever possible. At the end, a shape preserving approximation of a planar curve is discussed.

K~vwords." Approximation, Spline, Hyperbolic spline, Shape preserving spline.

1. Introduction

The problem of constructing a shape preserving approximation to given data has received quite an attention in the last years. Such approximations are particularly useful in C A D / C A M procedures. The general problem of determining a shape preserving approximation is nonlinear in its nature. Thus it has to be solved approximately, in one or more iterations. These iterations can be done in an interactive way by a human operator or by the computer itself. In this paper, we consider an algorithm that (approximately) does the job in one pass, and which is based upon the spline in tension. The idea of an interpolation spline under tension goes back to Schweikert [3]. He proposed it as a remedy in case the ordinary cubic interpolant does not preserve a shape of given data. His ideas were generalized by Sp~ith [4]. These generaliza-

7 Jernej Kozak was born in Ljubljana, Yugoslavia in 1946. He graduated in mathematics at the E.K. University of Ljuhljana in 1969. After graduation he became a full time programmer in the IMFM Computer Centre. He received a Master's degree in mathematics at the University of Zagreb (1972). In 1972 he was elected a lecturer at the department of mathematics of the E.K. University. He successfully defended tJ his PH.D. in 1977, and was elected Docent for Numerical Analysis and Computer Science at the E.K. University. In 1979/80 he was granted the Fulbright postdoctoral scholarship in Madison, Wisconsin, at the Mathematics Research Centre and University of Wisconsin. After his return he became more engaged in the cohaputer science, in particular in the theory of data structures and algorithms. Dr. Kozak has a lot of practical experience in practical numerical analysis and in programming in general. North-Holland Computers in Industry 7 (1986) 435-440

Fig. 1. Data connected by straight lines.

0166-3615/86/$3.50 ~2' 1986 Elsevier Science Publishers B.V. (North-Holland)

436

Short Note

(omputers in Industo"

J

[

/

\

/

/

/

/

Fig. 2. A first order interpolation.

Fig. 4. The taut spline interpolation.

tions allow local choices of the tension parameter, We consider this property an important one. Consider the example in Figs. 1-5. In Fig. l we see the data connected by straight lines; Fig. 2 gives an approximate contour obtained by approximation of the type:

In Fig. 3 the role of an approximating function is played by a cubic interpolatory spline, and in Fig. 4 de Boor's taut spline interpolation is applied [1]. The example clearly shows that cubic spline interpolation is quite superior to the case of Fig. 2, where all details are lost. However, around the nose and the chin the approximation is unsatisfactory. The local modification in taut spline recovers the shape of data also in this region but leaves the good parts unchanged. This is the reason why the local modifications as in de Boor's algorithm or in the procedure of this paper should be applied. The result of the algorithm of the paper is showed in Fig. 5. One of the disadvantages in the practical applications of splines in tension is the choice of the tension parameter(s). This paper can be considered as an elaboration on this remark.

E g ( ri)B,,4 .

J

2. The Shape Preserving Interpolatory Hyperbolic Spline Let ~':=(~)7-0 be a strictly increasing sequence of reals, with ~o := a, % := b. A piecewise exponential function f of order k with the breakpoint sequence ~- and the tension matrix ( p j,) R k'', restricted to [¢,, "ri~t), is defined by: Fig. 3. The cubic spline interpolation.

f[ ....... i ~ N ( L , ) ,

J, Kozak / Shape Preserving Approximation

Computers in Industry

437

and

/

p2sh ( p x ) ¢"(x) -

> o,

sh(p)

-

p

=

we conclude that f will have an inflection point in (0, 1) iff c3c 4 < 0. In our further discussion it will be easier to work with function and derivative values f, f ' at the endpoints rather then the coefficients q. But ¢(0) =¢(1) =0. Thus c,=f(0),

c 2 = / ( 1 ).

Further SO : = / " ( 0 )

~" - - C 1 ~1- C2 -- C3 -- 0/C4 '

s1:=f'(1) = -q Fig. 5. The hyperbolic spline interpolation.

+ c 2 + ac 3 + c a .

Since

p,,h(p)- sh(p) a:=a(p):=

sh(p)-p

1

1~ ]

where

p2

4[1

1~

~2,

p4

35 + ~ . + 57., + ' dxx

--

"'"

We shall restrict ourselves to piecewise exponential functions of order 4 [2], i.e. hyperbolic splines (of order 4):

and a2-1>3>0, we finally have

1 Li. . . . . dx 2 dx 2

pf

O
The reason for this is quite simple. A n increased n u m b e r of free parameters in one interval would disable us in keeping track of inflection points. Let us start with a simple case, n = 2, r 0 := O, r~ := 1. We choose the basis as in [4] and write a hyperbolic spline f in this special case as

f:=f(x)

:= Q n ( x ) + c2z/(1 - x ) "~-C31P(X) q- C4~b(1 -- X),

C4--

-- 1

s o + as 1

(oso+, ,~+1

+[0,1]/

)

.

Thus f will have no inflection point in (0, 1) iff

(~o + , ~ , ) ( ~ 0 + ~1) =<0 where the modified slopes si are defined by ~o:=so -[O,1]f,

~:=sl-[O,

1]f.

In these new terms f will have no inflection point if

with ~(x):=l-x,¢(x):=¢(x;

p)

,h ( e x ) - xsh ( p ) sh(p) -p ' where p, 0 __


f " = c3q/'( x ) + c44/'(1 - x ) ,

a+l-

-a+l'

with

~ 1 - - SO '

So far no conditions have been imposed on f yet. Suppose now that f agrees with a given g at

Short Note

438

(omputers in lndust O'

2 1~0 ]) is eliminated from the basis set. The interpolating hyperbolic spline will also continuously ~ 0 : = g ' ( 0 ) - [ 0 , a]g, ~ : = g ' ( 1 ) - [ 0 , l]f. be reduced to the straight line as g0 approaches 0. Thus the requirement that the only interpolating We do not want f to have an inflection point in function that goes through three colinear points is (0, 1) if there is no sign change in ([0, 0, 1]g, a straight line is satisfied. The argument for the [0, 1, a]g), i.e. other range of ~0 is obvious. ([0, 1]g - g'(O))( -[0, llg + g'(1)) = (-o~0)(: , ) > O. The Eq. 2 has a double root p = 0 for ~0=2. In this case the value of ¢0 satisfies Though one could try a more elaborate approach, the (monotone) Newton method will produce six 0
This conclusion holds for any p. We proceed to show how to choose p such that the inequalities (eq. 1) will not be possible. Note that: da

=

(sh(p)_ p) 2 ~-(-7.5+26-1)+

=p

(-9.7+28-1)+

p2

,.,2

p(,,,, )

TrY; 1/~<(

"

m = O, 1 .... Consider now the general case. The piecewise interpolatory hyperbolic spline f that agrees with a given g at r is on [r,, r , ~ ] given by

(ch(p)-l)(p+sh(p))-p2sh(p)

dp

( a( p{,,,~ )

p~,,, , l~ := p(,,,) _

... >0,

f:= t(x)

1 ( _ 2 . 6 + 25) + ~ ( ( _ 2.8 + 2v) + ..

~5

:=

0
('IiT~ AT,

as well as:

+£'3i+

l +O(p2e

-co.

gl --SO =

q,=g(r). C3/

At, O~i - 1

12'4,--

ai_

(2)

Any choice of ,b, p
!

[ 5, ffl -_~.i '

with

P).

Thus for any ~, ~ < ~o < 1, there exists a unique p such that a(p) a(p)+l

( .¥--T,

kT,

t

~(o) = 2, a(p)=p-

, Jr- C2i7 /

-- ' O/

12!4 ~ 0 .

The function that cannot bend fast enough at the endpoint where data show faster growth ( [sl I >

c2, = g(r,+ t),

s,+a,si+l _ [ r , r, ]g) a,+l

at, l

a,s,+s,+l a,+l

, , 1

-

+[7,. ~l]g')-

Here

,~,

:=

~( p, ),

where the value of Pi has to be determined so that f has no extraneous inflection point in [r,, r,~ ~]. The relevant quantities are now

.v0 : = ~ , - [~-,, ~,~:]g,

{3)

Vl := si~l - [r~, r,+l]g. If the given g is known twice at each of the data points r,, i.e. in the Hermit case. we have: •s'i = g'( r, ), all i, and the job is done. However, if the slopes sg are

Computers in Industry

J. Kozak / Shape Preserving Approximation

not known in advance but rather determined by

Then the equations

f e C(21[a, b],

f " ( r , - 0 ) = f " ( $ , + 0),

(4)

approximate values have to be used. Consider an approximation to ~0. If [r, ,, ~-,, ~-,+,]g=0,

become: 1

+

AT,fli I /~Tl[~i- 1 q._ - -

Ogi 1

a, 1 - 1

g:0 = 0, regardless of what the value of ~1 might be, since this will reduce f to the straight line in [r,, r,+l] as required in the case of Eq. 5. Thus only the data values g ( r , 1), g(r,), g(r,+l) can be used. From Eq. 3 we deduce that the only possible approximation to s, is of the form: 1"

i=1,2 ..... n-1

(5)

we want

2`[Ti

439

r, l g + (1 - X ) [ r , , r,+,]g.

The choice 2, = 1 can be found in [1]. Note that

fli l + f l ,

a,

1 ,8i l + f i ,

s,

+

Asfli 1 a,_l + l -~7-~+/~

a,_.,--1 [~, , ' ~ , ] g

AT, lfli a , + l + /~i_l+fli a7_--1 It,' 'r,+l]g' (8)

=1,2 ..... n-1.

kr, 2,'-kr,_l +At,

It is easy to verify that Eq. 8 is reduced to Eq. 9 [1], by the choice

gives an O(I,Ar I 2) approximation to g'(r,). This has to be compared with O( I Ar I) in the general case. Having this in mind we choose the approximate modified slopes ~0, Sl as:

p , = O,

at,

S0 := A,Tt_I-}- AT, ([Ti-I'

'Ti]g--

['Tt,

Ti+l]g ),

31-,

all i.

A remark has to be added. It is quite clear that the algorithm is not restricted only to interpolation problems. It can be applied provided that one is able to approximate the slope values at breakpoints as well as to estimate inflection points of a given function.

'~;I "-- aT,+, +aT, (['T'+I' T'+2]g--['T,, Tz+l]g ). 3. The Shape Preserving Approximation of a Curve

These two values are used to determine p,. We can now proceed to determine the spline. From Eq. 4 we obtain n - 1 equations, and there are several different ways to supply the two missing conditions. In particular, one might require

Let us consider a more general, hut important case: a shape preserving approximation of a curve. A straightforward approach along the lines of Section 2 gives us:

f "'(y 1-0)=f'"(r

X := X(I):= Cb'0 ~

1+0),

f'"(T,, 1 - - 0 ) = f ' " ( ' r , ,

]+0),

(6)

which is clearly an analog to a not-knot condition in the polynomial case. However, if the data slope values s 0 and s,, are not given (or approximated), one can discuss the shape of f only in [r 1, r, 1], and as a consequence in this case, pl=0,

p,, 1 = 0 .

fli "-

Pi sh ( P i )

ch(p,) - 1

, all i.

AT,

[t-ri

A~

)

; p' '

%
(7)

Though the equations, derived from Eq. 6 can be found in [4] in a more general context, we write them here for our particular case. Let us denote:

Jr- ('2,'0

>'=Y(t)::d,?[~

[t-%) C

+<,'0

t-r,

('r,+l - t

)

+ G, ¢ -57[-~, : p'

+d4i¢ T,+ I - t ] A.q ; P' ] "

~,

440

Short Note

Computers in lndusto

We shall throughout assume that t is a strictly increasing m o n o t o n e function of the natural parameter. Then an inflection point of r,

r:=r(t):=ly(t

[~',, ~+1] if there is no sign change in the data sequence

([~, ,, ~]x[~, 1, ~,, ,,+,ly,

. However, if the value of

is a zero of z,

. . . . .

z :=

SoxSII -- S1xSo~ ~

z(t)= ~

Recall again the simple case, n = 2, ~0 := 0, r~ := 1, with prescribed slopes now denoted by Sox, Sl,x, SOy, S l y ,

is small, then z(t) will have at most one zero in [0, 1]. In fact, if the parameter t is a natural one, we have

Since, on [0, 1],

z(t) = [0, 11~(

~ ( , ) = c,(~(t) + c~(~(1 - ,), as well as

aS°Y s+l " + ( 1 - - - t)') O~2 1 [ SOx + OtSlx

- [ 0 , l l y I -~F---~

))(t) = d 3 ~ ( t ) + d4~(1 - t), it follows that Eq. 10 is small if t is close to a natural parameter. Consequently, it is enough to keep track of inflection points in x and y separately. But then we can apply the procedure of the previous Section verbatim!

+(t)

ago~ + 21~

) /

+ (~ox~,-

(10)

{,? _ ~?)2 = ~2 + _ ~ 2

is given by JOy + o ~ l y

, a

d~i

d4

and further

Sox, gl~, Soy, st,,.

z(t)

C4

k2 +j,2 = 1,

and modified to

Then

c3

~,~0,,)(-+(t)(~O

- t)

+~(t)q)(a - t)).

(9)

References

However, the easy way now ends. The function -~(t)~(1

-t)

+ ~(t)q)(1 -t)

never vanishes in [0, 1]. As a consequence, z(t) can have two zeroes in [0, 1] what makes bookkeeping of inflection points in general a tedious task. Even so we could pursue the approach of Section 2, with an inflection being extraneous in

[1] C. de Boor (1978) A Practical Guide To Splines. Sprigerverlag, New York. [2] L. Schumaker (1981) Spline Functions: Basic Theory. John Wiley, New York. [3] D.G. Schweikert (1966) An Interpolating Curve Using a Spline in Tension. J. Math. Physics, 45. 312-317. [4] H. SpMh (1973) Spline-Algorithmen zur Konstruktion flatter Kuroen and Fliichen. R. Oldenbourg Verlag, Mfinchen.