- Email: [email protected]
On rational parametric curve approximation
Computer Aided North-Holland
COMAID
Geometric
Design
363
10 (1993) 363-377
342
On rational parametric curve approximation M.J. Pratt
*, R.J. Goult
** and L. Ye * * *
Dept. of Applied Computing Bedford MK43 OAL, UK
& Mathematics,
Presented in Oberwolfach, Received September 1992 Revised March 1993
June
Cranfield
Institute
of Technology,
Cranfield,
1992, by M.J. Pratt
Abstract Pratt, M.J., R.J. Goult and L. Ye, On rational Design 10 (1993) 363-377.
parametric
curve approximation,
Computer
Aided
Geometric
A new method is proposed for the rational approximation of functions. The error functional is chosen to avoid nonlinearity in the computation; the result is not a “best” approximation, but is usually good enough for many practical purposes. The same basic approach can be used for approximating either scalar or vector-valued functions. It is also possible to meet point and derivative constraints up to any desired order at curve endpoints. Keywords.
Rational
curves,
approximation;
offset curves;
geometric
continuity;
CAD data exchange.
1. Introduction Rational curves and surfaces have only become widely used in computer-aided design during the past ten years. However, the uses of rational functions were studied over many previous decades in the context of approximation theory, where they are generally found to provide better accuracy in approximation, for a given amount of computation, than pure polynomial functions [Ralston & Rabinowitz ‘851.Furthermore, they lend themselves well to the approximation of functions defined over infinite intervals or having singularities in the domain of interest [Hayes ‘701. The reasons for the use of rational curves and surfaces in CAGD are very different. The primary advantage claimed is that conic sections, quadric surfaces and surfaces of revolution having rational generators can all be represented exactly in terms of rational parametric geometry [Piegl & Tiller ‘87, Farin ‘891. Another potential advantage is the added freedom for shape design available through manipulation of the weights in the rational forms used [Faux & Pratt ‘79, Boehm ‘87, Farin ‘89, Su & Liu ‘891, though suitable methods are still being sought for presenting this freedom to the designer in such a way that it can be handled in an intuitive way. [Farin ‘911 contains a collection of recent research papers on rational curve and surface geometry. * Now at Design & Manufacturing Institute, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, * * Now at LMR Systems, 33 Filgrave, Newport Pagnell, Bucks MK16 9ET, UK. * * * Now at Capital Automation Inc., 6120 Harris Technology Boulevard, Charlotte, NC 28269, USA. 0167-8396/93/$06.00
0 1993 - Elsevier
Science
Publishers
B.V. All rights reserved
USA.
364
M.J. Pratt et al. / Rational approximation
The prevalence of rational representations of geometry in commercially available CAD systems has led to the provision of means for transfer of such geometry in the IGES and STEP standards for CAD data exchange [IGES ‘91, IS0 ‘881. The nature of such standards often influences decisions taken by developers of new systems, and for this reason we are likely to see more systems based on rational geometry in the future. Rational curves and surfaces also have disadvantages in CAGD, however. For example, problems arise due to poor parametrisation [Sederberg ‘861 and to lack of uniqueness in curve and surface representations [Farin ‘891. The space of rational splines is nonlinear, which often leads to a requirement for nonlinear calculations in the construction of curves or surfaces in constrained situations. Surfaces constructed as blends of smooth boundary curves may exhibit interior discontinuities of tangent plane. Computation of derivatives is much more expensive than for pure polynomial curve and surface representations, which makes for inefficiency in geometric interrogations such as the calculation of intersection curves. In view of these problems some authors have suggested that the current emphasis on rational geometry in CAGD reflects fashion rather than necessity [Farin ‘89, Dokken ‘911. On the whole the present authors agree with this sentiment. Nevertheless, it seems that rational geometry is here to stay, at least for a while. This fact has motivated the present study, which indicates that effectiveness in approximation for data exchange purposes may be another point in favour of rational representations in CAD. The problem of concern is the transfer of free-form curve and surface data between two CAD systems, specifically in the case when the receiving system implements rational geometry. Consider, for example, the Intergraph GMS system, which uses rational geometry of degree up to 15. Suppose we wish to transfer geometry into GMS from the Cisigraph system STRIM, which implements polynomial geometry up to degree 20. In general, then, it is not possible to effect an exact transfer. The STRIM curves and surfaces could be approximated by subdivision and the use of lower-degree polynomials, but here we examine the possibility, at least for curves, of generating rational rather than polynomial approximations. It will be shown that there are advantages in doing this, in that significantly better accuracy can be attained for a given complexity of approximation (in a sense to be defined in Section 3). Little has previously been published on rational approximation in a CAGD context. Hoschek [Hoschek ‘881 describes a method for rational cubic BCzier approximation to offset curves. Here the two inner control points and the weights for each curve segment are determined in terms of end constraints combined with a least-squares fit to computed points on the offset. At a late stage in the preparation of this paper the authors’ attention was drawn to unpublished work by Ambros [Ambros ‘891, who used a technique similar to that of this paper for the computation of rational cubic curve and surface approximations. Ambros computed the error functional in a different way, and his method, like that of Hoschek, is limited to approximants with cubic degree. By contrast, the method reported here is unrestricted as to degree, and is capable of matching derivative conditions of arbitrary order at the curve endpoints. It has application in the approximation of individual curve segments and also of their offsets. It may further be used in data reduction, i.e. the replacement of a piecewise defined curve by a single rational segment. Methods for computing best rational approximations ‘are in general nonlinear. Pade approximations provide an exception, but since these are based on Taylor expansions about the midpoint of the interval of interest the error does not have the uniform distribution appropriate for CAD purposes [Ralston & Rabinowitz ‘851. We here develop an efficient linear method yielding a good (but not best) rational L, approximation in a variety of situations. The result may be satisfactory as it stands for some purposes, or alternatively it will provide a good initial solution for the iterative computation of a best L, or L, rational approximation.
365
M.J. Pratt et al. / Rational approximation
2. L, function
approximation
using orthonormal
polynomials
Initially we consider the approximation of a scalar-valued function; the method to be described will later be generalised for use in vector-valued approximation. We will restrict attention to the interval [O, 11, and will work in terms of sets of polynomials Pi(t), i = 0, 1, 2,. . . ,where Pk(t) has degree k, orthonormal on this interval in the sense that
/0
‘W( t)P,( t)Pj( t)
dt = a,,,
(1)
where w(t) is a nonnegative weight function. For the moment, the weight function will be taken to be the unit function w(t) = 1. In this case the Pi(t) are the (normalised) Legendre polynomials, shifted and scaled from their usual interval of definition [ - 1, l] onto the interval [O, 11. It will be shown in Section 5 that alternative choices of weight functions permit approximation with positional and derivative constraints at the endpoints of the interval of interest. Initially, we briefly review the standard method of L, approximation in terms of the orthonormal polynomials Pi(t) as defined above. In approximating the function f(t) on [O, 11 we must first choose the degree m of the desired approximant P,(t), which is then written as a linear combination of the polynomials Pi(t) up to and including degree m: m
Pm(t)= The problem
c P&4(f).
(2)
i=o
is then to minimise
the error
functional
e(pm)= /ulW)-4mY dt =
jo’If2W - 2f(t)J’m(t)
with respect to the coefficients property (1) leads to
which may conveniently
+KiW} dt
pi. Substitution
be rewritten
(3)
of (2) into (3) and use of the orthonormality
as 2
e(f’,)
=i’.f’(t)
dt
+
I? 0
-
it
lof(t)Pi(f) (
(~f(tMf)
dt_Pi )
dt)‘.
(5)
0
Of the three terms, the first and third are independent of the pi. The second is nonnegative, and the error functional will clearly be minimised when this term is zero, for which the necessary conditions are pi =
lof(t)PiW dt,
i=o,
1 ,...,m.
(6)
366
M.J. Pratt et al. / Rational approximation
The minimum value of the error functional then follows immediately from (5) as
e(P,,J,,.=~lf2(f) dt-
EP~.
i=O
The approximant may be expressed in terms of basis functions other than the Pi(t), but then the corresponding procedure yields a set of linear equations, usually ill-conditioned, each equation in general containing all the unknown coefficients. The use of orthonormal polynomials avoids ill-conditioning by giving values for the coefficients directly with no requirement for solving equations. It also yields an explicit value for the minimum of e. Further, if a higher-degree approximation is desired there is no need for the entire calculation to be repeated; new higher-degree terms may be computed as required and appended to the original results. This useful property is known as permanence.
3. Approximation
by rational
functions
Suppose now that we wish to approximate rational function
the function f(t) not by a polynomial but by a
&n,(t) =f’rn(t)/Qn(t),
(8)
in which the numerator and denominator are polynomials of degree m, n respectively. It will be assumed that f(t) has no singularities on the interval [O, 11. The quantity m + n is known as the index of the rational function [Ralston & Rabinowitz ‘851, and it may be regarded as analogous to the degree of a pure polynomial. To take an example, Ro3(t), R,,(t), R,,(t) and R,,(t) all have index 3; the last of them is a pure polynomial. One coefficient of a rational approximation may always be chosen arbitrarily. This has the effect of scaling the numerator and denominator by the same factor. Hence there is the same number of unknown coefficients to be determined in any rational approximation of the same index (note, in this context, that Ro3(t), R,,(t) have respectively a constant numerator and a constant denominator). The computational effort of computing any rational approximation with a given index, since it requires the determination of the same number of coefficients, is roughly equivalent. All such approximations will be said to have the same complexity of approximation, which makes precise the statement on this topic occurring earlier in Section 1. If we try to use the method of Section 2 to compute the best L, approximation in terms of the unknown coefficients of Pm, Q, the problem proves to be a nonlinear one. Rational approximations are usually computed in practice by some iterative method, and success is most likely to be achieved if a good initial approximation is known. Remez’ Second Algorithm [Ralston & Rabinowitz ‘851 is one well-known iterative method, though it is used for computing best L,, rather than L, rational approximations. In this paper we avoid the inherent nonlinearity of the rational approximation problem by solving a related linear problem yielding a “good” result but not the best possible result. The chosen error functional is
E(fLn) = /‘{Q,WW) 0
-f’mW~2 dt.
The method therefore amounts to the simultaneous choice of a denominator such that Q,f “beh aves like a polynomial” and a numerator P,(t) which is a polynomial approximating it.
367
M.J. Pratt et al. / Rational approximation
Assume that we work in terms of the previously Pi(t), i = 0, 1, 2, . . . . We may then write pm(t)
=
Q,(f)
EPiDiCr)7
=
defined
set of orthonormal
basis functions
(10)
kqjBj(l).
0
0
Before embarking on the minimisation of E(R,,), recall that we can choose one of the coefficients arbitrarily. For present purposes we choose q,, = 1, which ensures that the denominator is not zero for t = 0; if it were, the approximant would have an undesirable singularity there. Following the procedure of Section 2, with the sole difference that f(t) is replaced by Q,(t)f(t), we find that 2
Then
the values
of the pi which minimise
pi=/“‘Q,(t)f(t)Pi(t) When
drj2-
E (~'Q,,(t)f(t)Pi(') 0
-
substituted
dt,
E(R,,), i=O,
(11)
for a given denominator
l,...,
m.
Q,(t),
are (12)
into (111, this gives
‘( ‘,,),i,
=
/o’e,“c
t)f2(
t>
dt
-
5
P?.
(13)
i=O
Equation (12) seems to suggest that the property of permanence mentioned in Section 2 still holds for the computation of the pi. In fact this is not so; as we shall see later, (21) implies that an increase in m requires Q,(t) to be recomputed. Then all the pi must be recalculated from scratch using (12). It remains to minimise E(R,,) with respect to the coefficients qj. TO this end we substitute the second of Eqs. (10) into (13) and equate to zero all its first partial derivatives is easily found from (121 with with respect to the qj, j = 1, 2,. . . , n. The derivative api/dqj substitution for Q,(t) in terms of the Pi(t) using (101. The result is 5
qk11f2(t)P,(t)Pk(‘)
k=O
dt-
j=l,2
EPi/‘f(t)Pi(t)Pj(t) i=o
o
,-..,
dt=O,
0
IZ.
(14)
We now define fij =
lOf(t)Pi(t )P,(r)
and note in passing that if f(t) other cases numerical quadrature forms for (12) and (14): ?
qjfrj-Pi=O,
dt,
gjk
=
I
‘.f2(t)Pj(t)Pk(
t,
dt,
(15)
0
is a polynomial the integrals may be evaluated exactly. In may be necessary. These definitions lead to more concise
i=O,
1,
j=l,
2 ,...,
. . . 7 m,
j=O (16) 2 k=O
qkgjk-
Epifij=O, i=O
n.
368
M.J. Pratt et al. / Rational approximation
This is a set of m + n + 1 linear equations in the unknown coefficients pO, pi,. . . , pm; since q. = 1 and hence the first term in each 9 2,. . . , qn. The system is not homogeneous, of the left-hand summations is a constant. For ease of manipulation the equations may be written in the following partitioned matrix form, where for future convenience the order of the two groups in (16) has first been reversed: 41,
[:
zT][:] =-[sf],
(17)
Here F is the (m + 1) X II matrix (fjj>, i = 0, 1,. . . , m; j = 1, 2,. . . , n, G is the II x n matrix C&j), i = 1, 2,. . .) n; j = 1, 2,. . .) ~1, Z is the (m + 1) X (m + 1) identity matrix and the superscript T denotes a transpose. The lowercase bold characters denote vectors, as follows: P= (PO, P,Y..,Pm)T~ f=(foo,
flOY.J-mO)T1
4= (s,, q2,..4JT, g=(g,,,
g*O,...,gn”)T.
(18)
Expansion of (17) restores in matrix form the two equations equivalent to (16), Gq=FTp-g
(19)
p=Fq+f.
(20)
and
Substitution of (20) into (19) leads to (G-FTF)q=FTf-g.
(21)
The elements of q may therefore be computed from (21) and then the elements of p from (20); this is computationally more efficient than a direct solution of the overall system (17). Note that it is the use of orthonormal basis functions which permits the overall solution to be found more efficiently by solving two smaller systems rather than a single system of order (m + IZ+ 1) X (m + n + 1). There is, as in Section 2, the added advantage that the calculation of the pi using (20) is direct and not ill-conditioned. However, p may not be computed until (21) has been solved for q, and the latter equation may prove to be singular. It is clear from (9) and (10) that this will occur whenever there is linear dependence amongst the set of m + n + 1 functions &,(t), Pi(t), . . . , P,(t), whose unknown coefficients are to be determined in conP,O>fO), &(t)./-(t), . . . , /3,(t)f(t) structing the approximation. In such a case the dependence implies that no unique solution exists, and that f(t) is itself a rational function with numerator and denominator having degree < m, q n respectively. Otherwise the stated set of functions will be linearly independent and a unique solution should be obtained. Since a meaningful solution may exist even in singular cases it is important to solve (21) by an appropriate method, and we have used singular value decomposition [Wilkinson ‘781. This allows solutions to be calculated for both nonsingular and singular examples. In singular cases any arbitrary coefficients qj arising are simply set to zero. Since we may obtain a singular system of equations in the circumstances described above, we may expect a system which is ill-conditioned (almost singular) when f(t) is a close approximation to a low-index rational function. Alternatively, ill-conditioning rather than singularity may arise even in the first case because of truncation errors arising from the use of numerical quadrature in the evaluation of the elements of the matrices G and F. Fortunately the singular value decomposition method handles ill-conditioning well, and even in the extreme case where numerical errors lead to a system which is actually inconsistent the
369
M.J. Pratt et al. / Rational approximation
method will give a best solution in the least-squares sense. No case of this last kind has yet been encountered in practice, however. One suggestion for detecting singular or ill-conditioned cases is to step systematically through sequences of increasing values of m and n. The occurrence of a very small residual for an index of approximation less than that originally intended will then indicate a possible degeneracy. For any solution of (17) the corresponding minimum value of the error functional is, from (13), E(R
4. Rational
mn ) mm = 2 0
approximation
iqjqkgjk0
EP,Z* 0
of a parametric
(22)
curve segment
Let r(t) = [x(t), y(t), z(t)], t E 10,11, be a parametric vector-valued approximation &z(t)
curve segment. We now seek a
=P,(t)/Q,(t)
(23)
which minimises the error functional E(&,)
= @(f)r(t)
-J’,(t)
I2 dt.
(24)
Using an obvious extension of the method of Section 3 we substitute Pm(t) =
F[Px,i~ Py,i,Pz,ilPiCt)7Qn( t) = IkqjP,(t) 3 0
(25)
0
where once again qt, = 1, and minimise the error functional with respect to the unknown coefficients. As previously it is convenient to express the values of the coefficients in terms of certain integrals involving the basis functions: Xij= / ‘x(t)Di(t)Dj(t) 0
dt>
zij=
dt,
Yij = /alY(t)Pi(t)Bj(t)
dt, (26)
/
‘z(t)Bi(t)Pj(t)
wjk =i1{x2(t)
+ty2(t)
+z’(t))Bj(t)p,(t)
dt.
Proceeding as previously, we once again obtain a system of linear equations, expressible in matrix form as
(27)
in this case
4
[;
--;
;
j,j
1 PX
PI’
PZ i
=-
(28)
370
Here
h4.J. Pratt et al. / Rational approximation
the new symbols
are defined
W=(W,~),
i=1,2
X=(xij),
i=O,
Y=(yij),
i=O,
Z=(zij),
i=O,l,...,
as follows:
,...,
n; j=l,2
,...,
n;
l,...,
m; j=l,2
,...,
n;
l,...,
m; j=l,Z
,...,
n;
m; j=l,2
,...,
n;
4= (41, 92,...4JT; Px = (Px,“, PxJ ,... &JT; Pz = (Pz,ot w=
P,,l....d,,,)T;
Pz,l,...7Pz,m)T;
(WlO,
Y = (Yoo,
P, = (P,,o,
wzo,...,
%ojT;
x=(x00,
T.
x,,,...,x,o)
7
T
Y1omYmlJT;
= =
(zoo,
The overall solution may once again be computed equations, obtained by expansion of (28) followed
=10,.
. *, z,o)
.
by solving smaller subsystems; by a little algebra, are
the relevant
(W-X-9-YTY-ZTZ)q=XTX+YTy+ZTz-w,
(29)
p,=xq
(30)
+x7
Py=Yq+Y,
(31)
p,=zq+z.
(32)
The results obtained specialise in an obvious way to planar parametric curves. Like Eq. (21), (29) may prove to be a singular system. Such cases are dealt with by the use of singular value decomposition as mentioned in the previous section. The derivation of the error bound in the parametric case follows that of equation (22); the result is E(R,,)min=
k 0
5. Rational
parametric
kqjqkwjkp 0
approximation
t(P,Z,i+P::l 0
with endpoint
+PZ,i).
(33)
constraints
Constrained approximation by vector-valued polynomials has been discussed in [Lachance ‘881, [Goult ‘891 and [Goult & Sherar ‘901. It is needed in CAD curve and surface data exchange to handle piecewise curves with specified derivative continuity between segments. The least-squares polynomial and rational approximations described in Sections 2, 3 and 4 have the disadvantage that the endpoints of the approximant do not coincide with those of the function being approximated. Then if these methods are used to approximate each segment of a smooth curve defined in a piecewise manner the composite approximant will not even possess Go continuity. The constrained methods cited above allow the approximation of piecewise defined curves (and also surfaces) by polynomials, with preservation of positional continuity and any desired degree of derivative continuity between segments (or patches). In this section the basis of Goult’s approach will be briefly reviewed. The method will then be extended to the generation of rational approximations, firstly with constrained end points but subsequently also with constrained end derivatives.
M.J. Pratt et al. / Rational approximation
371
Consider a set of polynomials orthonormal on [0, l] with respect to the weight function w(r) = t*(l - t>*. Denote this set of polynomials by r):(t), i = 0, 1, 2,. . . . Then the set of modified polynomials -~~+~(t)= t(1 - t>rl:(l) is clearly also orthonormal, but with respect to w(t) = 1; unlike the Pi(t), they all have the property that y(O) = ~(1) = 0. In fact the ri(t) are Jacobi polynomials [Abramowitz & Stegun ‘641, and the ri(t) are conveniently computed using an appropriate recurrence relation [Goult & Sherar ‘901. We now show how these constrained polynomials may be used in curve approximation with end constraints. Consider the scalar-valued function f(t), and define a new function 4(t) by 4(t) Note that endpoints Section 2 expressed
=f(t)
- {(I - t)f(O)
(34)
+rf@)J.
the term in curly braces is a linear interpolant between the values of f(t) at the of the interval of approximation, and that therefore 4(O) = 4(l) = 0. The method of is now used to approximate +(t> over [0, 11, with the modification that P,(t) is as a linear combination of the ri(t> rather than the P,(t):
(35) The summation starts with i = 2 since this is the lowest degree of basis function for which the endpoint constraints can be satisfied, as should be clear from the foregoing. Since all the ri(t) are zero at the endpoints the approximant agrees in value with 4(t) there. Then the function Pm(t) + (I - t)f(O> + 0Yl) is an approximation to the original function f(t) agreeing in value with it at the endpoints. Apart from the differences noted here the computation follows much the same lines as in Section 2. The same advantages result from the use of orthonormal polynomials. We now consider the extension of the method to rational approximation. The analysis of Section 3 needs little modification. As above, it is now $(t> rather than f(t) which is initially approximated, and P,(t) is expressed in terms of the ri(t> rather than the Pi(t). This is again appropriate since it is effectively Q,(t)+(t) which is being approximated by P,(t), and this, like 4(t) itself, is zero at the endpoints. There is however no apparent virtue in representing Q,(t) in terms of the ri(t), and in the approximation process we may therefore use two different bases: Pm(t> =
EPiYitt),
Q,(t)
=
?*jPj(t)*
(36)
0
2
On following through the analysis we arrive at a system of linear equations identical in form to (171,
[:
-fi’l[:] =-[;I>
(37)
though some of the symbols have modified definitions. Now F is the (m - 1) X n matrix (fij>, i=2,3 >. . *>m; j= 1,2 , . . .,y1, G is as before the n X n matrix (gij>, i = 1, 2,. . . , n; j = 1, 2,..., 12and Z is the (m - 1) x (m - 1) identity matrix. Vectors p and f have also changed, but q and g are as before: P’(P2,
P3,...,PmlT,
f= (fi0, f30Y.,f,0)T?
4=(41, g= (g,,,
q2Y..1qn)TY g20>...&O)T.
(38)
372
M.J. Pratt et al. / Rational approximation
The elements of f, F, g and G also have modified definitions: fij = i’+( t>ri( t)Pj( t) dt,
gij=IU’~‘(t)Pi(t)Pj(t)
dt.
(39)
The overall system of equations may as in Section 3 be written in terms of two subsystems, identical with (20) and (21). The error bound is given by (22) with p0 =p, = 0. We now turn to the constrained approximation of vector-valued functions. The method extends readily to this case, the analysis following that developed in Section 4. Once again we give the results for a 3D parametric curve. The linear system to be solved takes the same form as (28),
though there are some slight changes in the definitions of the elements of this equation: W=(wij),
i=l,2
,..., n;j=l,2
,..., n;
X=(X,~),
i=2,3
,...,
m;j=1,2
,..., n;
Y=(yij),
i=2,3
,...,
m; j=l,2
,..., n;
Z=(zij),
i=2,3
,...,
m; j=1,2
,..., n;
4=(41,
42,...,4JT; P, = (Py,2,
Px,3,...,Px,JT;
Px=(Px,2,
PY,3’ . . . . P,JTi
Pz=(Pr,z:P=;3,...,P;,,)T; w=
(W,“,
W20r...,
wno>‘;
x=(x,,,
X30,...,Xmo)T;
z=(z20,
z3o,...,z,o)
T
Y = (Yzo,
Y30~...,YmO)T~
.
The elements of w, W remain unchanged, but those of X, y, z, X, Y and 2 now involve the basis functions yi. If we define p(t)
= [P,(t),
p,(t),
p*(t)]
=r(t)
- (I - r)r(O) -@(I)1
(41)
then they are given by xij= / ‘P,(t)yi(t)Bj(t) 0 zlj =
/
olPz( t)ri(
dry
Yij = i’P,(t)Yi(t)Pj(t)
t>Pj(t) dt,
dt,
(42)
and
=k’{PZ(t) +Pc(t) +PI(t)}Pj(t)P/c(t) dt*
Wj/c
(43)
The overall system may be decomposed as in Section 4 in the form of Eqs. (29)-(321, for convenience in its solution. Once the rational approximant to p(t) is found the desired approximating curve is given by adding to it the linear interpolant represented by the last two
373
M.J. Pratt et al. / Rational approximation
terms in (41). The error bound is given by equation (33), slightly modified as for the scalar case. As previously, the results specialise easily to the approximation of planar curves. In this and the previous example in this section the systems of equations (37) and (40) decompose into simpler systems by virtue of the representation of P,(t), or P,(t), in terms of the constrained orthonormal basis functions ri(t). The basis functions Pi(t) used in the representation of Q,(t) play no role in the simplification of the linear system. Evaluation of the integrals in (39), (42) and (43) will in fact be simplified if the Pi(t) are replaced in the analysis by the simple monomial basis I’, though the effects of this replacement on the stability and accuracy of the computations have not so far been investigated.
6. Rational
parametric
approximation
with endpoint
and derivative
constraints
The final development in this paper is the rational approximation of vector-valued functions subject to both positional and derivative constraints at the endpoints. Suppose that we desire the approximant to have the same tangent vector as the original function r(t) at these points. A natural approach is to define a new set of polynomials A,(t), orthonormal with respect to the weight function w(t) = t4(1 - t>4 (these are further cases of Jacobi polynomials). Then the modified polynomials Ai( t) = t2( 1 - t)‘A,(t)
(44)
will be orthonorrnal with respect to w(t) = 1 and will have their values and first derivatives zero at the endpoints of the interval [O, 11. The function corresponding to p(t) in (41) must be obtained by subtracting from r(t) a function sharing the same endpoints and end tangent vectors. The simplest choice is now to employ a cubic interpolant having the desired properties, namely the well-known Ferguson/Hermite cubic curve segment [Faux & Pratt ‘791. The resulting difference has null position and tangent vectors for t = 0, 1, and it is therefore appropriate to approximate it as a sum of the h(t), since this gives the same properties. The details are otherwise much the same as before. The method described above may be generalised in an obvious manner to enable the matching of derivative vectors up to any desired order k at the endpoints of the interval of approximation. The most straightforward choice for the function initially subtracted is the polynomial of degree 2k + 1 having matching end conditions. If this is denoted by J2k+l(t), the overall approximant is PVZ(t)/Q,(t)
+J2k+l
(f)
=
kW
+Q,(t>J,,+,(t>}/Q,(t),
(45)
sothat the final rational approximant has its numerator of degree max(m, n + 2k + 1). Since P,(r) is constructed in terms of constrained basis functions as described above, we must have m a 2k. Further, if we wish to prescribe the degree of the numerator on the right-hand side of (45) as being m, then we must also require the stronger condition m>2k+n+l,
(46) will be larger than originally
otherwise the index of the overall rational approximation specified. It is worthy of note that the derivatives of the approximant at its end-points agree in both direction and magnitude with those of the original curve function. Then if a composite curve is approximated segment by segment the continuity characteristics of that curve will be reproduced exactly at the junction points provided that k in the above analysis is given an appropriate value. In practice, it will not often be necessary to choose k > 2; k = 2 will give either C2 or G2 continuity, consistent with the nature of the original piecewise curve.
M.J. Pratt et al. / Rational approximation
374
Fig. 1.
7. Examples
The examples are concerned with the approximation f(t)
= (1t
of the scalar function
25t2)-1’2
on the interval [ - 1, l¬e that the function is symmetrical about t = 0. Initially we compute two unconstrained approximations to f(t), the first a polynomial of degree 6 and the second a rational function having degree 3 in both numerator and denominator. Both approximating functions have index 6, and their complexity is therefore the same (i.e. we have seven free parameters to calculate in each case). It is clear from Fig. 1 that the rational approximation has a maximum error less than half that of the polynomial approximation, and that as we move away from t = 0 the rational form gives greatly superior results. The rational parametric approximation resulting from setting x = t, y =f(t) gives virtually indistinguishable results and is therefore not plotted separately. In Fig. 2 we show the results of constraining the approximation to have the same end-points as f(t) at t = & 1. This requires the degree of the numerator to be raised. For
Fig. 2.
375
M.J. Pratt et al. / Rational approximation
*c----__
I--_
I
I’
-----’
,‘-
=.__/
,’
-\
‘\ ‘.___--
.. /_----__
,---. --___*
I\
Fig. 3.
comparability with the previous results we consider the case with seven coefficients to be determined; the degree of the denominator is again 3, but that of the numerator is now 5. The result is clearly better than that of the unconstrained rational approximation in Fig. 1, though it is found that in some cases the increased degree of the numerator contributes less towards overall accuracy than to meeting the end constraints. Figure 3 shows the unconstrained approximation of f(t) by a rational function with degree 5 in both numerator and denominator. The errors are now such that the two curves are hardly distinguishable. Superimposed is the polynomial approximation of degree 8, which still exhibits marked oscillations about the curve of .f(t>.
8. Discussion
and conclusions
A method has been proposed which is capable of giving good rational approximations to either scalar or vector-valued functions without the necessity for solving nonlinear equations. Remez’ algorithm, referred to earlier, determines best L, rational approximations through the solution of the same number of equations, one for each unknown coefficient, but in this case the equations are nonlinear and must be solved iteratively. Our tradeoff is therefore increased computational efficiency against optimal&y of approximation. There are situations where the former is more important, and our experience has been that the rational approximations given by the present, linear, method are significantly better in many cases than polynomial approximations requiring the same order of computational effort. The method is currently at an early stage of refinement, and several aspects of it require further study. These include the following. 1. Expression of a bound on the true error functional e(R,,) in terms of the modified error functional E(R,,) as defined in Sections 3 and 4. 2. Investigation of methods for constraining Q,(t) to be nonnegative on [O, l]-or, better still, to restrict its deviation from the value unity. One suggestion for achieving this is to use the basis functions ri(t> defined in Section 5, and to represent Q,(t) by
Qn(t) = 1 + CqjYj(‘), 2
376
M.J. Pratt et al. / Rational approximation
which applies point constraints at the endpoints of the interval of approximation. This has not yet been tried, but it seems likely that such constraints on the denominator will require a correspondingly higher degree in the numerator to achieve the same accuracy of approximation. 3. Elucidation of the effect of different choices of basis functions on the stabilty of the solution for the coefficients sj. 4. In the context of the constrained approximation methods of Sections 5 and 6, study of the effect of subtracting different interpolants satisfying the same end conditions in constructing the subsidiary function which is initially approximated. 5. Clearly, for CAGD purposes, the method needs to be extended for the rational approximation of surfaces. We have described extensions to the basic method which permit the generation of approximations satisfying point constraints and derivative constraints of arbitrary order at the endpoints of the interval of interest. In a CAGD context, this allows the segment-by-segment approximation of piecewise defined curves in such a way that the continuity characteristics of the original curve are retained by the overall approximant. The method may also be used to approximate a sequence of original curve segments by a single rational curve. It also lends itself to the approximation of offset curves and any other kinds of procedurally defined curves which permit points and derivatives to be evaluated for arbitrary parameter values. In approximating parametric curves the error functional is based on distance between points on the original and approximating curves which have the same parameter values. The approximation is therefore not based on purely geometric principles, as it would be, for instance, if the criterion were based on the Hausdorff distance between the two curves. On the other hand the criterion used leads to more efficient procedures than would then be possible, and in many cases of the transfer of engineering data between CAD systems there is a strong argument for approximating not only the geometry but also its parametrisation. The method introduced here achieves this effect. Support for the virtues of a linear approach to rational approximation is given by the previously mentioned work [Ambros ‘891. His investigations were limited to approximants with cubic degree in both numerator and denominator, and he compared a linear method with two alternative nonlinear methods. As might be expected the linear approach was much the most efficient computationally, by a factor ranging from 5 to 20 in the curve approximation case. However, although the accuracy attained was usually comparable with that of the nonlinear methods based on the true error functional (cf. Eq. (3) but with a rational approximant), some examples were found where the error was much larger. This reinforces the need for further enhancement of the method as suggested above. Some initial work has been done on the extension of the technique to the construction of rational curves providing an approximate fit to point data. This is based on the principle of orthogonality under summation rather than, as here, integration. The results will be reported more fully in a future paper. 9. Acknowledgement Lin Ye gratefully acknowledges the financial support of the British Council during his visit to Cranfield Institute of Technology. 10. References Abramowitz, M. and Stegun, LA., eds., (19651, Handbook of Mathematical Functions, Dover, New York. Ambros, M. (1989), Approximation von BLzierfltichen durch rationale BCzierfllchen vom Grad (3,3), Master’s Fachbereich Mathematik, Technische Hochschule Darmstadt, Germany.
thesis,
M.J. Pratt et al. / Rational approximation
377
Boehm, W. (1987), Rational geometric splines, Computer Aided Geometric Design 4, 67-78. Dokken, T., Ytrehus, A.M. and Skytt, V. (1991) The role of NURBS in geometric modelling and CAD/CAM, in: Krause, F.-L. and Jansen, H. eds., Advanced Geometric Mode&g for Engineering Applications (Proc. IFIP WG52/GI Internat. Symp., Berlin, Nov. 1989), North-Holland, Amsterdam. Farin, G. (1989) Rational curves and surfaces, in: Lyche, T. and Schumaker L.L., eds., Mathematical Methods in Computer Aided Geometric Design, Academic Press, New York. Farin, G., ed., (1991) NUMBS for Curve and Surface Design, SIAM, Philadelphia, PA. Faux, I.D. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood, Chichester, UK. Goult, R.J. (1989) Parametric curve and surface approximation, in: Handscomb, D.C., ed., The Mathematics of Surfaces III (Proc. IMA Conf., Oxford, Sept. 1988), Oxford Univ. Press, Oxford. Goult, R.J. and Sherar, P.A., eds., (1990), Improving the Performance of Neutral File Data Transfers, ESPRIT Research Reports Series (Project 322 - CAD Interfaces, Vol. 6). Springer, Berlin. Hayes, J.G. (1970) Numerical Approximations to Functions and Data, Athlone Press. Hoschek, J. (1988) Spline approximation of offset curves, Computer Aided Geometric Design 5, 33-40. IGES (1991) Initial Graphics Exchange Specification (IGES), Version 5.1, Sept. 1991, NCGA Technical Services & Standards, Fairfax, VA, USA. IS0 (1993), DIS 10303-42: Product Data Representation and Exchange, Part 42 - Integrated Generic Resources: Geometrical and Topological Representation, Draft International Standard, International Standards Organization. Lachance, M.A. (1988) Chebyshev economisation for parametric surfaces, Computer Aided Geometric Design 5, 195-208. Piegl, L. and Tiller, W. (1987) Curve and surface constructions using rational B-splines, Computer Aided Design 19, 499-502. Ralston, A. and Rabinowitz, P. (1985) A First Course in Numerical Analysis, McGraw-Hill, New York. Sederberg, T.W. (1986), Improperly parametrized rational curves, Computer Aided Geometric Design 3, 67-75. Su, B. and Liu, D. (1989) Computational Geometry-Curve and Surface Modelling, Academic Press, New York. Wilkinson, J.H. (1978), Singular-value decomposition-basic aspects, in: Jacobs, D.A.H., ed., Numerical SoftwareNeeds and Availability, Academic Press, New York.
COMAID
Geometric
Design
363
10 (1993) 363-377
342
On rational parametric curve approximation M.J. Pratt
*, R.J. Goult
** and L. Ye * * *
Dept. of Applied Computing Bedford MK43 OAL, UK
& Mathematics,
Presented in Oberwolfach, Received September 1992 Revised March 1993
June
Cranfield
Institute
of Technology,
Cranfield,
1992, by M.J. Pratt
Abstract Pratt, M.J., R.J. Goult and L. Ye, On rational Design 10 (1993) 363-377.
parametric
curve approximation,
Computer
Aided
Geometric
A new method is proposed for the rational approximation of functions. The error functional is chosen to avoid nonlinearity in the computation; the result is not a “best” approximation, but is usually good enough for many practical purposes. The same basic approach can be used for approximating either scalar or vector-valued functions. It is also possible to meet point and derivative constraints up to any desired order at curve endpoints. Keywords.
Rational
curves,
approximation;
offset curves;
geometric
continuity;
CAD data exchange.
1. Introduction Rational curves and surfaces have only become widely used in computer-aided design during the past ten years. However, the uses of rational functions were studied over many previous decades in the context of approximation theory, where they are generally found to provide better accuracy in approximation, for a given amount of computation, than pure polynomial functions [Ralston & Rabinowitz ‘851.Furthermore, they lend themselves well to the approximation of functions defined over infinite intervals or having singularities in the domain of interest [Hayes ‘701. The reasons for the use of rational curves and surfaces in CAGD are very different. The primary advantage claimed is that conic sections, quadric surfaces and surfaces of revolution having rational generators can all be represented exactly in terms of rational parametric geometry [Piegl & Tiller ‘87, Farin ‘891. Another potential advantage is the added freedom for shape design available through manipulation of the weights in the rational forms used [Faux & Pratt ‘79, Boehm ‘87, Farin ‘89, Su & Liu ‘891, though suitable methods are still being sought for presenting this freedom to the designer in such a way that it can be handled in an intuitive way. [Farin ‘911 contains a collection of recent research papers on rational curve and surface geometry. * Now at Design & Manufacturing Institute, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, * * Now at LMR Systems, 33 Filgrave, Newport Pagnell, Bucks MK16 9ET, UK. * * * Now at Capital Automation Inc., 6120 Harris Technology Boulevard, Charlotte, NC 28269, USA. 0167-8396/93/$06.00
0 1993 - Elsevier
Science
Publishers
B.V. All rights reserved
USA.
364
M.J. Pratt et al. / Rational approximation
The prevalence of rational representations of geometry in commercially available CAD systems has led to the provision of means for transfer of such geometry in the IGES and STEP standards for CAD data exchange [IGES ‘91, IS0 ‘881. The nature of such standards often influences decisions taken by developers of new systems, and for this reason we are likely to see more systems based on rational geometry in the future. Rational curves and surfaces also have disadvantages in CAGD, however. For example, problems arise due to poor parametrisation [Sederberg ‘861 and to lack of uniqueness in curve and surface representations [Farin ‘891. The space of rational splines is nonlinear, which often leads to a requirement for nonlinear calculations in the construction of curves or surfaces in constrained situations. Surfaces constructed as blends of smooth boundary curves may exhibit interior discontinuities of tangent plane. Computation of derivatives is much more expensive than for pure polynomial curve and surface representations, which makes for inefficiency in geometric interrogations such as the calculation of intersection curves. In view of these problems some authors have suggested that the current emphasis on rational geometry in CAGD reflects fashion rather than necessity [Farin ‘89, Dokken ‘911. On the whole the present authors agree with this sentiment. Nevertheless, it seems that rational geometry is here to stay, at least for a while. This fact has motivated the present study, which indicates that effectiveness in approximation for data exchange purposes may be another point in favour of rational representations in CAD. The problem of concern is the transfer of free-form curve and surface data between two CAD systems, specifically in the case when the receiving system implements rational geometry. Consider, for example, the Intergraph GMS system, which uses rational geometry of degree up to 15. Suppose we wish to transfer geometry into GMS from the Cisigraph system STRIM, which implements polynomial geometry up to degree 20. In general, then, it is not possible to effect an exact transfer. The STRIM curves and surfaces could be approximated by subdivision and the use of lower-degree polynomials, but here we examine the possibility, at least for curves, of generating rational rather than polynomial approximations. It will be shown that there are advantages in doing this, in that significantly better accuracy can be attained for a given complexity of approximation (in a sense to be defined in Section 3). Little has previously been published on rational approximation in a CAGD context. Hoschek [Hoschek ‘881 describes a method for rational cubic BCzier approximation to offset curves. Here the two inner control points and the weights for each curve segment are determined in terms of end constraints combined with a least-squares fit to computed points on the offset. At a late stage in the preparation of this paper the authors’ attention was drawn to unpublished work by Ambros [Ambros ‘891, who used a technique similar to that of this paper for the computation of rational cubic curve and surface approximations. Ambros computed the error functional in a different way, and his method, like that of Hoschek, is limited to approximants with cubic degree. By contrast, the method reported here is unrestricted as to degree, and is capable of matching derivative conditions of arbitrary order at the curve endpoints. It has application in the approximation of individual curve segments and also of their offsets. It may further be used in data reduction, i.e. the replacement of a piecewise defined curve by a single rational segment. Methods for computing best rational approximations ‘are in general nonlinear. Pade approximations provide an exception, but since these are based on Taylor expansions about the midpoint of the interval of interest the error does not have the uniform distribution appropriate for CAD purposes [Ralston & Rabinowitz ‘851. We here develop an efficient linear method yielding a good (but not best) rational L, approximation in a variety of situations. The result may be satisfactory as it stands for some purposes, or alternatively it will provide a good initial solution for the iterative computation of a best L, or L, rational approximation.
365
M.J. Pratt et al. / Rational approximation
2. L, function
approximation
using orthonormal
polynomials
Initially we consider the approximation of a scalar-valued function; the method to be described will later be generalised for use in vector-valued approximation. We will restrict attention to the interval [O, 11, and will work in terms of sets of polynomials Pi(t), i = 0, 1, 2,. . . ,where Pk(t) has degree k, orthonormal on this interval in the sense that
/0
‘W( t)P,( t)Pj( t)
dt = a,,,
(1)
where w(t) is a nonnegative weight function. For the moment, the weight function will be taken to be the unit function w(t) = 1. In this case the Pi(t) are the (normalised) Legendre polynomials, shifted and scaled from their usual interval of definition [ - 1, l] onto the interval [O, 11. It will be shown in Section 5 that alternative choices of weight functions permit approximation with positional and derivative constraints at the endpoints of the interval of interest. Initially, we briefly review the standard method of L, approximation in terms of the orthonormal polynomials Pi(t) as defined above. In approximating the function f(t) on [O, 11 we must first choose the degree m of the desired approximant P,(t), which is then written as a linear combination of the polynomials Pi(t) up to and including degree m: m
Pm(t)= The problem
c P&4(f).
(2)
i=o
is then to minimise
the error
functional
e(pm)= /ulW)-4mY dt =
jo’If2W - 2f(t)J’m(t)
with respect to the coefficients property (1) leads to
which may conveniently
+KiW} dt
pi. Substitution
be rewritten
(3)
of (2) into (3) and use of the orthonormality
as 2
e(f’,)
=i’.f’(t)
dt
+
I? 0
-
it
lof(t)Pi(f) (
(~f(tMf)
dt_Pi )
dt)‘.
(5)
0
Of the three terms, the first and third are independent of the pi. The second is nonnegative, and the error functional will clearly be minimised when this term is zero, for which the necessary conditions are pi =
lof(t)PiW dt,
i=o,
1 ,...,m.
(6)
366
M.J. Pratt et al. / Rational approximation
The minimum value of the error functional then follows immediately from (5) as
e(P,,J,,.=~lf2(f) dt-
EP~.
i=O
The approximant may be expressed in terms of basis functions other than the Pi(t), but then the corresponding procedure yields a set of linear equations, usually ill-conditioned, each equation in general containing all the unknown coefficients. The use of orthonormal polynomials avoids ill-conditioning by giving values for the coefficients directly with no requirement for solving equations. It also yields an explicit value for the minimum of e
3. Approximation
by rational
functions
Suppose now that we wish to approximate rational function
the function f(t) not by a polynomial but by a
&n,(t) =f’rn(t)/Qn(t),
(8)
in which the numerator and denominator are polynomials of degree m, n respectively. It will be assumed that f(t) has no singularities on the interval [O, 11. The quantity m + n is known as the index of the rational function [Ralston & Rabinowitz ‘851, and it may be regarded as analogous to the degree of a pure polynomial. To take an example, Ro3(t), R,,(t), R,,(t) and R,,(t) all have index 3; the last of them is a pure polynomial. One coefficient of a rational approximation may always be chosen arbitrarily. This has the effect of scaling the numerator and denominator by the same factor. Hence there is the same number of unknown coefficients to be determined in any rational approximation of the same index (note, in this context, that Ro3(t), R,,(t) have respectively a constant numerator and a constant denominator). The computational effort of computing any rational approximation with a given index, since it requires the determination of the same number of coefficients, is roughly equivalent. All such approximations will be said to have the same complexity of approximation, which makes precise the statement on this topic occurring earlier in Section 1. If we try to use the method of Section 2 to compute the best L, approximation in terms of the unknown coefficients of Pm, Q, the problem proves to be a nonlinear one. Rational approximations are usually computed in practice by some iterative method, and success is most likely to be achieved if a good initial approximation is known. Remez’ Second Algorithm [Ralston & Rabinowitz ‘851 is one well-known iterative method, though it is used for computing best L,, rather than L, rational approximations. In this paper we avoid the inherent nonlinearity of the rational approximation problem by solving a related linear problem yielding a “good” result but not the best possible result. The chosen error functional is
E(fLn) = /‘{Q,WW) 0
-f’mW~2 dt.
The method therefore amounts to the simultaneous choice of a denominator such that Q,
367
M.J. Pratt et al. / Rational approximation
Assume that we work in terms of the previously Pi(t), i = 0, 1, 2, . . . . We may then write pm(t)
=
Q,(f)
EPiDiCr)7
=
defined
set of orthonormal
basis functions
(10)
kqjBj(l).
0
0
Before embarking on the minimisation of E(R,,), recall that we can choose one of the coefficients arbitrarily. For present purposes we choose q,, = 1, which ensures that the denominator is not zero for t = 0; if it were, the approximant would have an undesirable singularity there. Following the procedure of Section 2, with the sole difference that f(t) is replaced by Q,(t)f(t), we find that 2
Then
the values
of the pi which minimise
pi=/“‘Q,(t)f(t)Pi(t) When
drj2-
E (~'Q,,(t)f(t)Pi(') 0
-
substituted
dt,
E(R,,), i=O,
(11)
for a given denominator
l,...,
m.
Q,(t),
are (12)
into (111, this gives
‘( ‘,,),i,
=
/o’e,“c
t)f2(
t>
dt
-
5
P?.
(13)
i=O
Equation (12) seems to suggest that the property of permanence mentioned in Section 2 still holds for the computation of the pi. In fact this is not so; as we shall see later, (21) implies that an increase in m requires Q,(t) to be recomputed. Then all the pi must be recalculated from scratch using (12). It remains to minimise E(R,,) with respect to the coefficients qj. TO this end we substitute the second of Eqs. (10) into (13) and equate to zero all its first partial derivatives is easily found from (121 with with respect to the qj, j = 1, 2,. . . , n. The derivative api/dqj substitution for Q,(t) in terms of the Pi(t) using (101. The result is 5
qk11f2(t)P,(t)Pk(‘)
k=O
dt-
j=l,2
EPi/‘f(t)Pi(t)Pj(t) i=o
o
,-..,
dt=O,
0
IZ.
(14)
We now define fij =
lOf(t)Pi(t )P,(r)
and note in passing that if f(t) other cases numerical quadrature forms for (12) and (14): ?
qjfrj-Pi=O,
dt,
gjk
=
I
‘.f2(t)Pj(t)Pk(
t,
dt,
(15)
0
is a polynomial the integrals may be evaluated exactly. In may be necessary. These definitions lead to more concise
i=O,
1,
j=l,
2 ,...,
. . . 7 m,
j=O (16) 2 k=O
qkgjk-
Epifij=O, i=O
n.
368
M.J. Pratt et al. / Rational approximation
This is a set of m + n + 1 linear equations in the unknown coefficients pO, pi,. . . , pm; since q. = 1 and hence the first term in each 9 2,. . . , qn. The system is not homogeneous, of the left-hand summations is a constant. For ease of manipulation the equations may be written in the following partitioned matrix form, where for future convenience the order of the two groups in (16) has first been reversed: 41,
[:
zT][:] =-[sf],
(17)
Here F is the (m + 1) X II matrix (fjj>, i = 0, 1,. . . , m; j = 1, 2,. . . , n, G is the II x n matrix C&j), i = 1, 2,. . .) n; j = 1, 2,. . .) ~1, Z is the (m + 1) X (m + 1) identity matrix and the superscript T denotes a transpose. The lowercase bold characters denote vectors, as follows: P= (PO, P,Y..,Pm)T~ f=(foo,
flOY.J-mO)T1
4= (s,, q2,..4JT, g=(g,,,
g*O,...,gn”)T.
(18)
Expansion of (17) restores in matrix form the two equations equivalent to (16), Gq=FTp-g
(19)
p=Fq+f.
(20)
and
Substitution of (20) into (19) leads to (G-FTF)q=FTf-g.
(21)
The elements of q may therefore be computed from (21) and then the elements of p from (20); this is computationally more efficient than a direct solution of the overall system (17). Note that it is the use of orthonormal basis functions which permits the overall solution to be found more efficiently by solving two smaller systems rather than a single system of order (m + IZ+ 1) X (m + n + 1). There is, as in Section 2, the added advantage that the calculation of the pi using (20) is direct and not ill-conditioned. However, p may not be computed until (21) has been solved for q, and the latter equation may prove to be singular. It is clear from (9) and (10) that this will occur whenever there is linear dependence amongst the set of m + n + 1 functions &,(t), Pi(t), . . . , P,(t), whose unknown coefficients are to be determined in conP,O>fO), &(t)./-(t), . . . , /3,(t)f(t) structing the approximation. In such a case the dependence implies that no unique solution exists, and that f(t) is itself a rational function with numerator and denominator having degree < m, q n respectively. Otherwise the stated set of functions will be linearly independent and a unique solution should be obtained. Since a meaningful solution may exist even in singular cases it is important to solve (21) by an appropriate method, and we have used singular value decomposition [Wilkinson ‘781. This allows solutions to be calculated for both nonsingular and singular examples. In singular cases any arbitrary coefficients qj arising are simply set to zero. Since we may obtain a singular system of equations in the circumstances described above, we may expect a system which is ill-conditioned (almost singular) when f(t) is a close approximation to a low-index rational function. Alternatively, ill-conditioning rather than singularity may arise even in the first case because of truncation errors arising from the use of numerical quadrature in the evaluation of the elements of the matrices G and F. Fortunately the singular value decomposition method handles ill-conditioning well, and even in the extreme case where numerical errors lead to a system which is actually inconsistent the
369
M.J. Pratt et al. / Rational approximation
method will give a best solution in the least-squares sense. No case of this last kind has yet been encountered in practice, however. One suggestion for detecting singular or ill-conditioned cases is to step systematically through sequences of increasing values of m and n. The occurrence of a very small residual for an index of approximation less than that originally intended will then indicate a possible degeneracy. For any solution of (17) the corresponding minimum value of the error functional is, from (13), E(R
4. Rational
mn ) mm = 2 0
approximation
iqjqkgjk0
EP,Z* 0
of a parametric
(22)
curve segment
Let r(t) = [x(t), y(t), z(t)], t E 10,11, be a parametric vector-valued approximation &z(t)
curve segment. We now seek a
=P,(t)/Q,(t)
(23)
which minimises the error functional E(&,)
= @(f)r(t)
-J’,(t)
I2 dt.
(24)
Using an obvious extension of the method of Section 3 we substitute Pm(t) =
F[Px,i~ Py,i,Pz,ilPiCt)7Qn( t) = IkqjP,(t) 3 0
(25)
0
where once again qt, = 1, and minimise the error functional with respect to the unknown coefficients. As previously it is convenient to express the values of the coefficients in terms of certain integrals involving the basis functions: Xij= / ‘x(t)Di(t)Dj(t) 0
dt>
zij=
dt,
Yij = /alY(t)Pi(t)Bj(t)
dt, (26)
/
‘z(t)Bi(t)Pj(t)
wjk =i1{x2(t)
+ty2(t)
+z’(t))Bj(t)p,(t)
dt.
Proceeding as previously, we once again obtain a system of linear equations, expressible in matrix form as
(27)
in this case
4
[;
--;
;
j,j
1 PX
PI’
PZ i
=-
(28)
370
Here
h4.J. Pratt et al. / Rational approximation
the new symbols
are defined
W=(W,~),
i=1,2
X=(xij),
i=O,
Y=(yij),
i=O,
Z=(zij),
i=O,l,...,
as follows:
,...,
n; j=l,2
,...,
n;
l,...,
m; j=l,2
,...,
n;
l,...,
m; j=l,Z
,...,
n;
m; j=l,2
,...,
n;
4= (41, 92,...4JT; Px = (Px,“, PxJ ,... &JT; Pz = (Pz,ot w=
P,,l....d,,,)T;
Pz,l,...7Pz,m)T;
(WlO,
Y = (Yoo,
P, = (P,,o,
wzo,...,
%ojT;
x=(x00,
T.
x,,,...,x,o)
7
T
Y1omYmlJT;
= =
(zoo,
The overall solution may once again be computed equations, obtained by expansion of (28) followed
=10,.
. *, z,o)
.
by solving smaller subsystems; by a little algebra, are
the relevant
(W-X-9-YTY-ZTZ)q=XTX+YTy+ZTz-w,
(29)
p,=xq
(30)
+x7
Py=Yq+Y,
(31)
p,=zq+z.
(32)
The results obtained specialise in an obvious way to planar parametric curves. Like Eq. (21), (29) may prove to be a singular system. Such cases are dealt with by the use of singular value decomposition as mentioned in the previous section. The derivation of the error bound in the parametric case follows that of equation (22); the result is E(R,,)min=
k 0
5. Rational
parametric
kqjqkwjkp 0
approximation
t(P,Z,i+P::l 0
with endpoint
+PZ,i).
(33)
constraints
Constrained approximation by vector-valued polynomials has been discussed in [Lachance ‘881, [Goult ‘891 and [Goult & Sherar ‘901. It is needed in CAD curve and surface data exchange to handle piecewise curves with specified derivative continuity between segments. The least-squares polynomial and rational approximations described in Sections 2, 3 and 4 have the disadvantage that the endpoints of the approximant do not coincide with those of the function being approximated. Then if these methods are used to approximate each segment of a smooth curve defined in a piecewise manner the composite approximant will not even possess Go continuity. The constrained methods cited above allow the approximation of piecewise defined curves (and also surfaces) by polynomials, with preservation of positional continuity and any desired degree of derivative continuity between segments (or patches). In this section the basis of Goult’s approach will be briefly reviewed. The method will then be extended to the generation of rational approximations, firstly with constrained end points but subsequently also with constrained end derivatives.
M.J. Pratt et al. / Rational approximation
371
Consider a set of polynomials orthonormal on [0, l] with respect to the weight function w(r) = t*(l - t>*. Denote this set of polynomials by r):(t), i = 0, 1, 2,. . . . Then the set of modified polynomials -~~+~(t)= t(1 - t>rl:(l) is clearly also orthonormal, but with respect to w(t) = 1; unlike the Pi(t), they all have the property that y(O) = ~(1) = 0. In fact the ri(t) are Jacobi polynomials [Abramowitz & Stegun ‘641, and the ri(t) are conveniently computed using an appropriate recurrence relation [Goult & Sherar ‘901. We now show how these constrained polynomials may be used in curve approximation with end constraints. Consider the scalar-valued function f(t), and define a new function 4(t) by 4(t) Note that endpoints Section 2 expressed
=f(t)
- {(I - t)f(O)
(34)
+rf@)J.
the term in curly braces is a linear interpolant between the values of f(t) at the of the interval of approximation, and that therefore 4(O) = 4(l) = 0. The method of is now used to approximate +(t> over [0, 11, with the modification that P,(t) is as a linear combination of the ri(t> rather than the P,(t):
(35) The summation starts with i = 2 since this is the lowest degree of basis function for which the endpoint constraints can be satisfied, as should be clear from the foregoing. Since all the ri(t) are zero at the endpoints the approximant agrees in value with 4(t) there. Then the function Pm(t) + (I - t)f(O> + 0Yl) is an approximation to the original function f(t) agreeing in value with it at the endpoints. Apart from the differences noted here the computation follows much the same lines as in Section 2. The same advantages result from the use of orthonormal polynomials. We now consider the extension of the method to rational approximation. The analysis of Section 3 needs little modification. As above, it is now $(t> rather than f(t) which is initially approximated, and P,(t) is expressed in terms of the ri(t> rather than the Pi(t). This is again appropriate since it is effectively Q,(t)+(t) which is being approximated by P,(t), and this, like 4(t) itself, is zero at the endpoints. There is however no apparent virtue in representing Q,(t) in terms of the ri(t), and in the approximation process we may therefore use two different bases: Pm(t> =
EPiYitt),
Q,(t)
=
?*jPj(t)*
(36)
0
2
On following through the analysis we arrive at a system of linear equations identical in form to (171,
[:
-fi’l[:] =-[;I>
(37)
though some of the symbols have modified definitions. Now F is the (m - 1) X n matrix (fij>, i=2,3 >. . *>m; j= 1,2 , . . .,y1, G is as before the n X n matrix (gij>, i = 1, 2,. . . , n; j = 1, 2,..., 12and Z is the (m - 1) x (m - 1) identity matrix. Vectors p and f have also changed, but q and g are as before: P’(P2,
P3,...,PmlT,
f= (fi0, f30Y.,f,0)T?
4=(41, g= (g,,,
q2Y..1qn)TY g20>...&O)T.
(38)
372
M.J. Pratt et al. / Rational approximation
The elements of f, F, g and G also have modified definitions: fij = i’+( t>ri( t)Pj( t) dt,
gij=IU’~‘(t)Pi(t)Pj(t)
dt.
(39)
The overall system of equations may as in Section 3 be written in terms of two subsystems, identical with (20) and (21). The error bound is given by (22) with p0 =p, = 0. We now turn to the constrained approximation of vector-valued functions. The method extends readily to this case, the analysis following that developed in Section 4. Once again we give the results for a 3D parametric curve. The linear system to be solved takes the same form as (28),
though there are some slight changes in the definitions of the elements of this equation: W=(wij),
i=l,2
,..., n;j=l,2
,..., n;
X=(X,~),
i=2,3
,...,
m;j=1,2
,..., n;
Y=(yij),
i=2,3
,...,
m; j=l,2
,..., n;
Z=(zij),
i=2,3
,...,
m; j=1,2
,..., n;
4=(41,
42,...,4JT; P, = (Py,2,
Px,3,...,Px,JT;
Px=(Px,2,
PY,3’ . . . . P,JTi
Pz=(Pr,z:P=;3,...,P;,,)T; w=
(W,“,
W20r...,
wno>‘;
x=(x,,,
X30,...,Xmo)T;
z=(z20,
z3o,...,z,o)
T
Y = (Yzo,
Y30~...,YmO)T~
.
The elements of w, W remain unchanged, but those of X, y, z, X, Y and 2 now involve the basis functions yi. If we define p(t)
= [P,(t),
p,(t),
p*(t)]
=r(t)
- (I - r)r(O) -@(I)1
(41)
then they are given by xij= / ‘P,(t)yi(t)Bj(t) 0 zlj =
/
olPz( t)ri(
dry
Yij = i’P,(t)Yi(t)Pj(t)
t>Pj(t) dt,
dt,
(42)
and
=k’{PZ(t) +Pc(t) +PI(t)}Pj(t)P/c(t) dt*
Wj/c
(43)
The overall system may be decomposed as in Section 4 in the form of Eqs. (29)-(321, for convenience in its solution. Once the rational approximant to p(t) is found the desired approximating curve is given by adding to it the linear interpolant represented by the last two
373
M.J. Pratt et al. / Rational approximation
terms in (41). The error bound is given by equation (33), slightly modified as for the scalar case. As previously, the results specialise easily to the approximation of planar curves. In this and the previous example in this section the systems of equations (37) and (40) decompose into simpler systems by virtue of the representation of P,(t), or P,(t), in terms of the constrained orthonormal basis functions ri(t). The basis functions Pi(t) used in the representation of Q,(t) play no role in the simplification of the linear system. Evaluation of the integrals in (39), (42) and (43) will in fact be simplified if the Pi(t) are replaced in the analysis by the simple monomial basis I’, though the effects of this replacement on the stability and accuracy of the computations have not so far been investigated.
6. Rational
parametric
approximation
with endpoint
and derivative
constraints
The final development in this paper is the rational approximation of vector-valued functions subject to both positional and derivative constraints at the endpoints. Suppose that we desire the approximant to have the same tangent vector as the original function r(t) at these points. A natural approach is to define a new set of polynomials A,(t), orthonormal with respect to the weight function w(t) = t4(1 - t>4 (these are further cases of Jacobi polynomials). Then the modified polynomials Ai( t) = t2( 1 - t)‘A,(t)
(44)
will be orthonorrnal with respect to w(t) = 1 and will have their values and first derivatives zero at the endpoints of the interval [O, 11. The function corresponding to p(t) in (41) must be obtained by subtracting from r(t) a function sharing the same endpoints and end tangent vectors. The simplest choice is now to employ a cubic interpolant having the desired properties, namely the well-known Ferguson/Hermite cubic curve segment [Faux & Pratt ‘791. The resulting difference has null position and tangent vectors for t = 0, 1, and it is therefore appropriate to approximate it as a sum of the h(t), since this gives the same properties. The details are otherwise much the same as before. The method described above may be generalised in an obvious manner to enable the matching of derivative vectors up to any desired order k at the endpoints of the interval of approximation. The most straightforward choice for the function initially subtracted is the polynomial of degree 2k + 1 having matching end conditions. If this is denoted by J2k+l(t), the overall approximant is PVZ(t)/Q,(t)
+J2k+l
(f)
=
kW
+Q,(t>J,,+,(t>}/Q,(t),
(45)
sothat the final rational approximant has its numerator of degree max(m, n + 2k + 1). Since P,(r) is constructed in terms of constrained basis functions as described above, we must have m a 2k. Further, if we wish to prescribe the degree of the numerator on the right-hand side of (45) as being m, then we must also require the stronger condition m>2k+n+l,
(46) will be larger than originally
otherwise the index of the overall rational approximation specified. It is worthy of note that the derivatives of the approximant at its end-points agree in both direction and magnitude with those of the original curve function. Then if a composite curve is approximated segment by segment the continuity characteristics of that curve will be reproduced exactly at the junction points provided that k in the above analysis is given an appropriate value. In practice, it will not often be necessary to choose k > 2; k = 2 will give either C2 or G2 continuity, consistent with the nature of the original piecewise curve.
M.J. Pratt et al. / Rational approximation
374
Fig. 1.
7. Examples
The examples are concerned with the approximation f(t)
= (1t
of the scalar function
25t2)-1’2
on the interval [ - 1, l¬e that the function is symmetrical about t = 0. Initially we compute two unconstrained approximations to f(t), the first a polynomial of degree 6 and the second a rational function having degree 3 in both numerator and denominator. Both approximating functions have index 6, and their complexity is therefore the same (i.e. we have seven free parameters to calculate in each case). It is clear from Fig. 1 that the rational approximation has a maximum error less than half that of the polynomial approximation, and that as we move away from t = 0 the rational form gives greatly superior results. The rational parametric approximation resulting from setting x = t, y =f(t) gives virtually indistinguishable results and is therefore not plotted separately. In Fig. 2 we show the results of constraining the approximation to have the same end-points as f(t) at t = & 1. This requires the degree of the numerator to be raised. For
Fig. 2.
375
M.J. Pratt et al. / Rational approximation
*c----__
I--_
I
I’
-----’
,‘-
=.__/
,’
-\
‘\ ‘.___--
.. /_----__
,---. --___*
I\
Fig. 3.
comparability with the previous results we consider the case with seven coefficients to be determined; the degree of the denominator is again 3, but that of the numerator is now 5. The result is clearly better than that of the unconstrained rational approximation in Fig. 1, though it is found that in some cases the increased degree of the numerator contributes less towards overall accuracy than to meeting the end constraints. Figure 3 shows the unconstrained approximation of f(t) by a rational function with degree 5 in both numerator and denominator. The errors are now such that the two curves are hardly distinguishable. Superimposed is the polynomial approximation of degree 8, which still exhibits marked oscillations about the curve of .f(t>.
8. Discussion
and conclusions
A method has been proposed which is capable of giving good rational approximations to either scalar or vector-valued functions without the necessity for solving nonlinear equations. Remez’ algorithm, referred to earlier, determines best L, rational approximations through the solution of the same number of equations, one for each unknown coefficient, but in this case the equations are nonlinear and must be solved iteratively. Our tradeoff is therefore increased computational efficiency against optimal&y of approximation. There are situations where the former is more important, and our experience has been that the rational approximations given by the present, linear, method are significantly better in many cases than polynomial approximations requiring the same order of computational effort. The method is currently at an early stage of refinement, and several aspects of it require further study. These include the following. 1. Expression of a bound on the true error functional e(R,,) in terms of the modified error functional E(R,,) as defined in Sections 3 and 4. 2. Investigation of methods for constraining Q,(t) to be nonnegative on [O, l]-or, better still, to restrict its deviation from the value unity. One suggestion for achieving this is to use the basis functions ri(t> defined in Section 5, and to represent Q,(t) by
Qn(t) = 1 + CqjYj(‘), 2
376
M.J. Pratt et al. / Rational approximation
which applies point constraints at the endpoints of the interval of approximation. This has not yet been tried, but it seems likely that such constraints on the denominator will require a correspondingly higher degree in the numerator to achieve the same accuracy of approximation. 3. Elucidation of the effect of different choices of basis functions on the stabilty of the solution for the coefficients sj. 4. In the context of the constrained approximation methods of Sections 5 and 6, study of the effect of subtracting different interpolants satisfying the same end conditions in constructing the subsidiary function which is initially approximated. 5. Clearly, for CAGD purposes, the method needs to be extended for the rational approximation of surfaces. We have described extensions to the basic method which permit the generation of approximations satisfying point constraints and derivative constraints of arbitrary order at the endpoints of the interval of interest. In a CAGD context, this allows the segment-by-segment approximation of piecewise defined curves in such a way that the continuity characteristics of the original curve are retained by the overall approximant. The method may also be used to approximate a sequence of original curve segments by a single rational curve. It also lends itself to the approximation of offset curves and any other kinds of procedurally defined curves which permit points and derivatives to be evaluated for arbitrary parameter values. In approximating parametric curves the error functional is based on distance between points on the original and approximating curves which have the same parameter values. The approximation is therefore not based on purely geometric principles, as it would be, for instance, if the criterion were based on the Hausdorff distance between the two curves. On the other hand the criterion used leads to more efficient procedures than would then be possible, and in many cases of the transfer of engineering data between CAD systems there is a strong argument for approximating not only the geometry but also its parametrisation. The method introduced here achieves this effect. Support for the virtues of a linear approach to rational approximation is given by the previously mentioned work [Ambros ‘891. His investigations were limited to approximants with cubic degree in both numerator and denominator, and he compared a linear method with two alternative nonlinear methods. As might be expected the linear approach was much the most efficient computationally, by a factor ranging from 5 to 20 in the curve approximation case. However, although the accuracy attained was usually comparable with that of the nonlinear methods based on the true error functional (cf. Eq. (3) but with a rational approximant), some examples were found where the error was much larger. This reinforces the need for further enhancement of the method as suggested above. Some initial work has been done on the extension of the technique to the construction of rational curves providing an approximate fit to point data. This is based on the principle of orthogonality under summation rather than, as here, integration. The results will be reported more fully in a future paper. 9. Acknowledgement Lin Ye gratefully acknowledges the financial support of the British Council during his visit to Cranfield Institute of Technology. 10. References Abramowitz, M. and Stegun, LA., eds., (19651, Handbook of Mathematical Functions, Dover, New York. Ambros, M. (1989), Approximation von BLzierfltichen durch rationale BCzierfllchen vom Grad (3,3), Master’s Fachbereich Mathematik, Technische Hochschule Darmstadt, Germany.
thesis,
M.J. Pratt et al. / Rational approximation
377
Boehm, W. (1987), Rational geometric splines, Computer Aided Geometric Design 4, 67-78. Dokken, T., Ytrehus, A.M. and Skytt, V. (1991) The role of NURBS in geometric modelling and CAD/CAM, in: Krause, F.-L. and Jansen, H. eds., Advanced Geometric Mode&g for Engineering Applications (Proc. IFIP WG52/GI Internat. Symp., Berlin, Nov. 1989), North-Holland, Amsterdam. Farin, G. (1989) Rational curves and surfaces, in: Lyche, T. and Schumaker L.L., eds., Mathematical Methods in Computer Aided Geometric Design, Academic Press, New York. Farin, G., ed., (1991) NUMBS for Curve and Surface Design, SIAM, Philadelphia, PA. Faux, I.D. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood, Chichester, UK. Goult, R.J. (1989) Parametric curve and surface approximation, in: Handscomb, D.C., ed., The Mathematics of Surfaces III (Proc. IMA Conf., Oxford, Sept. 1988), Oxford Univ. Press, Oxford. Goult, R.J. and Sherar, P.A., eds., (1990), Improving the Performance of Neutral File Data Transfers, ESPRIT Research Reports Series (Project 322 - CAD Interfaces, Vol. 6). Springer, Berlin. Hayes, J.G. (1970) Numerical Approximations to Functions and Data, Athlone Press. Hoschek, J. (1988) Spline approximation of offset curves, Computer Aided Geometric Design 5, 33-40. IGES (1991) Initial Graphics Exchange Specification (IGES), Version 5.1, Sept. 1991, NCGA Technical Services & Standards, Fairfax, VA, USA. IS0 (1993), DIS 10303-42: Product Data Representation and Exchange, Part 42 - Integrated Generic Resources: Geometrical and Topological Representation, Draft International Standard, International Standards Organization. Lachance, M.A. (1988) Chebyshev economisation for parametric surfaces, Computer Aided Geometric Design 5, 195-208. Piegl, L. and Tiller, W. (1987) Curve and surface constructions using rational B-splines, Computer Aided Design 19, 499-502. Ralston, A. and Rabinowitz, P. (1985) A First Course in Numerical Analysis, McGraw-Hill, New York. Sederberg, T.W. (1986), Improperly parametrized rational curves, Computer Aided Geometric Design 3, 67-75. Su, B. and Liu, D. (1989) Computational Geometry-Curve and Surface Modelling, Academic Press, New York. Wilkinson, J.H. (1978), Singular-value decomposition-basic aspects, in: Jacobs, D.A.H., ed., Numerical SoftwareNeeds and Availability, Academic Press, New York.