- Email: [email protected]
Parametric approximation of data using ODR splines
COMPUTER AIDED GEOMETRIC DESIGN
Computer Aided Geometric Design 11 (1994) 247-267
Parametric
approximation Samuel
of data using ODR splines
P. Marin a,*, Philip
W. Smith b
a Mathematics Department, GM Research and Development Center, Warren, MI 48090-9055, USA b 14311 Kellywood, Houston, TX 77079, USA
Received June 1990; revised August 1991
Abstract We consider curve (that is, and numerical arises naturally
the problem of approximating data in space that arises from a geometric a curve with no explicit parameterization given). We treat theoretical aspects of a parametric orthogonal distance regression problem which in the context of geometric data fitting.
Key words: Orthogonal B-spline curves
distance
regression;
Parametric
approximation;
Least squares;
1. Introduction Given m ordered data points in Rd, {A}, i = 1,. . . , m, we want to find a curve which approximates the data and has a relatively simple form. By a curve we mean a function g which maps the unit interval continuously into Rd. Consider, for example, the data shown in Fig. 1. The points here were measured using a coordinate measuring machine along a cross section of a spiral surface of an air compressor component. The goal is to fit a parametric curve to the data and thus analytically represent the shape of the cross section. We will measure the approximation power of curves fitted to such data by computing E(g,f):=min 1
eJ/f;-g(Ti)J/‘:~iE it
[O,l],i=
l,...,m}.
(1)
I
We emphasize that this is not the classical measure of discrepancy which would require us to choose a priori the points where the errors are to be measured. Of course, if no restriction is placed on g, then we could just interpolate the *Corresponding author. Email: [email protected]. 0167-8396/94/$07.00 @ 1994 Elsevier Science B.V. All rights reserved 0167-8396(93)E0029-D
SSDI
248
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
D
Y (mm) 0.0 0
D
-50.0 _j
1
-50.0
0.0
50.0
x (mm) Fig. 1. Sample CMM data for parametric
curve in R2.
data at ri = (i - 1 )/ (m - 1). This would be unacceptable in most situations where the data are noisy. Furthermore, it would mean that the function g would become more complex as m grows. We avoid this issue of complexity by restricting the coordinate functionals of g to lie in a finite dimensional spline subspace. This leads us to the following fundamental approximation problem which we consider in this paper: inf{E(s,f):Pi(s)ES~,i=
l,...,
d}.
(2)
Here, SF is the subspace of kth order spline functions with knot sequence t (see de Boor (1978)) and Pi is the projection operator onto the ith coordinate. Thus s is a curve in Rd with each coordinate functional a spline in Sf. In the following we assume that both k and t are known and fixed. We also require that the interval spanned by the knot sequence t contains [0,1 ] so that the spline functions are defined over the target parameterization interval. We call solutions to (2) (if they exist) ODR splines. The ODR stands for orthogonal distance regression. This is discussed in (Boggs et al., 1987) where an algorithm is presented and analyzed to solve the ODR problem for explicitly represented curves and surfaces. The parametric problem posed and studied in our paper represents a variant which is not directly addressed by the methods outlined in (Boggs et al., 1987). In (Boggs et al., 1989), an extension of the methods in (Boggs et al., 1987) to multiple response data is described. This provides a capability to handle parametric problems. However, the objective function used in (Boggs et al., 1989 ) differs from (2) and therefore leads to a different geometric interpretation of the results. We say more about
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
-50.0 _! -50.0
249
0.0 x (mm)
Fig. 2. Ordinary least squares fit to data from Fig. 1.
this in Section 2.3. The parametric problem we consider has been previously posed by Plass and Stone ( 1983) and by Hoschek ( 1988) and references given there. These references present numerical approaches together with some examples from CAD applications. Additional examples and comparisons to other approaches are given by O’Dell ( 1985 ). We provide additional insights pertaining to the existence, uniqueness, characterization, and the approximation power of the solutions to these problems. A core element in the parametric ODR problem is that of data parameterization. For related literature covering other aspects of data parameterization for parametric spline curves we refer to (Foley and Nielson, 1989; Lee, 1989; Marin, 1984; Pinkus, 1988; Rademacher and Scherer, 1990; Scherer and Smith, 1989; Scherer, 1987; Toepfer, 198 1). There is a clear difference between the approximation problem above and the usual method for solving the parametric curve fitting problem. The standard approach to the data fitting problem requires the selection of a subspace as above in which the coordinate functionals must lie. Additionally, points 0 6 rr d *.. d rm Q 1 are selected. Now the standard least squares problem is solved (3) This standard least squares approach may be applied to the data of Fig. 1 by first choosing, for example, {7i, i = 1, . . . , m} to be normalized accumulated chord length and SF to be piecewise cubits defined on [0, 1 ] with seven
250
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
Y (mm)
0.0 -
-50.0 J -50.0
0.0
50.0
x (mm) Fig. 3. ODR spline fit to data from Fig. 1.
uniformly spaced interior breakpoints. The results are shown in Fig. 2. This standard least squares problem is linear. The price of linearity, however, is the inflexibility of prescribing the points ri. One of the principal advantages of the ODR spline problem (2) is that it incorporates the selection of {r} as part of its solution and thus alleviates this inflexibility. To see what might be gained by resorting to an ODR spline fit, we present, in Fig. 3, the parametric, piecewise cubic solution of (2) for the data shown in Fig. 1. The spline space SF is the same here as in the ordinary least squares case discussed above. These results, however, show a noticeable improvement over those displayed in Fig. 2. In the remainder of the paper we expand our discussion of this approach to parametric curve fitting according to the following outline. In Section 2 we describe a solution algorithm and present several numerical examples which provide additional illustrations of the capabilities and limitations of ODR splines. Section 3 contains theoretical results addressing questions of existence, uniqueness, and the characterization of solutions of these problems. In Section 4 we comment on underlying approximation issues. We conclude with a brief summary in Section 5.
2. Numerical examples The main purpose of this section is to present several examples which illustrate both the pitfalls and the power of ODR splines versus the classical
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
251
least squares method. Solutions for the examples chosen are obtained using a simple iterative technique. The method, which for completeness we describe below, is essentially that presented in (Hoschek, 1988). At the end of this section we comment on the connection between the data fitting problem as we have described it and orthogonal distance regression. 2.1. Algorithm description The starting point is the observation that if s* (q) solves (2), then it is necessary that s* also solve the ordinary least squares problem (3) with ri, i = l,..., m chosen according to
s*(7i)ll
Il_h-
= ,@, 1l.h- s*(rl)lI. . .
(4)
The algorithm described in (Hoschek, 1988) proceeds by iterating between problems (3) and (4). Starting with an initial parameterization {r(O)} and a spline s(O), (s(O) solves (3) with ri = ri”, i = 1,. . .,m), s(I) are generated which satisfy, respectively,
iterates
{rj”}
and
and
min
{
2 IIf;: - scr)(711,)l12 : Pj(s(l))
E
s:).
(6)
i=l
In the implementation discussed in (Hoschek, 1988)) the univariate minimization problem (5 ) is solved approximately by first constructing, at q = $-I) a local linear approximation to the parametric curve s (I-l ) (v ). With ,11-l,‘. m (5) replaced by its linear approximate, an updated parameter value, ,!I) from I ) is determined
7% I
(S+lqp),fi =
p
+
v
[
- su-l)(p))
~I,cr-w(7!~-9~~2 I
1.
Here, v is a scalar search variable, initially equal to 1, which is halved until the predicted rj” satisfies [IA - s”-‘)(s~“)ll < llfi - ~(~-‘)(r~‘-“)ll. In most cases we used this approach. However, in a few demanding examples where the fitted curve developed regions of high curvature, this local method was not appropriate. To proceed in these situations, we relied on a simple global, uniform search to solve (5 ) . Normally, we expect that the fitted curve will preserve the order of the data implied by the initial parametrization. That is, if r;, i = 1,. . . , m is the final parametrization determined by solving (2)) then we desire that
252
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
It is clear that the minimization problem (5) may not yield solutions r!l’, i = 1, . . . , m which satisfy the monotonicity condition (7). The likeliness of this happening depends, in part, on the choice of the spline subspace. For example, if we chose to solve the ODR problem for the data of Fig. 1 using parametric splines whose projections are linear polynomials, then solutions of (5) would not form a monotone sequence as desired. The remedy here is to increase the order k and refine the knot sequence t so that the spline subspace 5’: can more closely approximate the data. Our limited experience suggests that, for large enough k and fine enough t, the monotonicity condition (7) will be satistied as the iteration proceeds. We can, however, envision pathological examples where this will not be the case. In some cases it may be desirable to include the constraints (7) explicitly in the problem formulation. It may also be possible to solve such order constrained ODR spline problems using an iterative technique similar to that described above for the problem without order constraints. To do this we would change the update procedure for the vector r (0 of parameter values. Instead of solving the sequence of univariate minimization problems given in (5 ) to obtain r!” ,I = l,..., m, these values would be derived as the solution of the multivariate problem
The iteration would proceed by cycling between (8 ) and (6 ) . Although we have not carried out numerical experiments to test this approach, we do derive, in Section 3.4, the necessary conditions for a solution in this setting. Without a guarantee of existence of solutions for this parametric ODR problem (see Section 3), it is difficult to comment extensively on the convergence of the numerical algorithms. One useful point can be made, however. The iterates {r(“} and {s(l)} constructed according to (5 ) and (6 ) can be shown to satisfy E(s(‘-‘),f)
2 E(s’“,f),
1 = l)...
(9)
E (s, f ) is defined by ( 1). Thus the sequence of objective values {E (s(l), f )} is monotone decreasing, bounded below and hence converges to a value E*. As we will illustrate later, the spline coefficients associated with
where
the iterates {.s(‘)} may, however, grow without bound. We now consider several examples to highlight the capabilities limitations of ODR splines.
and the
2.2. Examples In the examples that follow, we seek ODR spline fits to four two-dimensional problems. In each case, the data is initially parameterized on the unit interval [0,1 ] and the underlying spline space $ is chosen so that the interval spanned by its knots strictly contains [0, 11. The only constraints imposed in solving the
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
0.5
I/
0.5
0.0
1.0
253
0.5
0.0
1.0
7
7
Fig. 4. Ordinary least squares fit to circle data.
Fig. 5. ODR spline fit to circle data.
minimization problem (5) are that 0 d ~j” d 1, i = 1,. . . , m. The graphs that accompany these examples show the parametric ODR spline, s*(q), plotted for q E [r;,z&].
Example 1. For our first example we construct the solution of problem (2) when the data {J;:}, i = 1,. . . , m, consist of 30 points in [w2uniformly spaced on the unit circle. The spline space consists of cubic B-splines with 11 knots t1 = t2 = t3 = t4 = -0.01, t5 = 0.245,
t(j = 0.5,
t8 = t9 = tlo = tll
= 1.01.
t7 = 0.755,
254
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
-ll.Ol
.
-10.0
. 1 -9.0
,
,
, , -8.6
,
,
( -7.0
109 (Cco)
Fig. 6. Piecewise linear interpolant of sample reaction rate data.
The data is parameterized (0)
T1
=
0,
,y
=
initially by normalized
T!O)+ r-1
IN - “LIII Ci”=,llfi
-jj_lI1’
accumulated
chord length,
.
E= 2,***,m*
The ordinary least squares fit to this data is given in Fig. 4. A plot of the fitted parametric spline is shown along with a plot of curvature. The same information is given in Fig. 5, but for the ODR spline fit to the data. Notable differences between Figs. 4 and 5 concern the behavior of the parametric curves near the data point ( 1,O) (no periodicity condition was imposed in either case) and the magnitude of the observed oscillations in the curvature plots. These are significantly lower in the case of the ODR spline. We note that periodicity can be imposed in the spline subspace without changing the solution methodology. We chose not to do so in this example to illustrate how well the base method captures geometric properties of the data.
Example 2. The second example
illustrates a potential application of these methods to laboratory data fitting. The situation we have in mind concerns a reacting chemical system which may possess multiple steady states. To obtain data for the example, we have sampled one of the model curves given in (Herz and Marin, 1980). The abscissa, Cco, of each data point is an input concentration of carbon monoxide in a catalytic system. The ordinate, R, is the steady state oxidation rate achieved by the system. In Fig. 6 we have plotted (log (Cc0 ),log (R ) ) at each sampled point and, to establish order, the piecewise linear interpolant of the data. Since multiple reaction rates are attained at a single input concentration, parametric curve fitting methods must be used.
255
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
-8.0
0
-8.0
OS0
(b)
(a) -9.c
-9.0
lo9 w
lo9 (W
-10.0
!
-11.o
-10.0
-8.0
-9.0
-7.0
-10.0
-9.0
-8.0
-7.0
log (Cco)
b2(Cco)
-8.0
-9.0
log (R) -10.0
-11.0 I_. -10.0
-9.0
-8.0
-7.0
~og~cco)
Fig. 7. Ordinary
least squares fits to reaction rate data. (a) One interior knot. (b) Two interior knots. (c) Three interior knots.
With the data parameterized by normalized accumulated chord length we show, in Fig. 7 (a), the ordinary least squares fit using parametric cubic splines having breakpoints at -0.01, 0.50, and 1.Ol. With only one interior breakpoint, the data are not well fit. One way to improve the fit is to introduce more breakpoints. In Figs. 7(b) and 7(c) we show the ordinary least squares fit obtained using, respectively, two and three uniformly spaced interior breakpoints in the interval (-0.01, 1.01). The curve in Fig. 7 (b) is still a rather poor fit to the data. It can be improved somewhat by more carefully positioning the two interior breakpoints. The curve shown in Fig. 7 (c) (three interior knots) is better, but still fails to model the behavior near the peak in the data. The ODR spline fitting technique offers an alternative procedure that is not so sensitive to the number or location of breakpoints. In fact, we can very effectively fit the reaction rate data with a parametric cubic spline curve
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
256
-8.0
-8.0
1
-9.0
-9.0.
‘09 CR)
‘og CR)
-10.0.
-10.0.
-1l.OI
-11.0I_
-10.0
1
-9.0
-8.0
-7.0
-10.0
109 (Cc,)
-9.0
-8.0
-7.0
~og(cco)
-8.0
-11.0A -10.0
-9.0
-8.0
-7.0
'09Gzo)
Fig. 8. ODR fit to reaction
rate data (one interior knot). (c) 1000 iterations.
(a) 10 iterations.
(b) 100 iterations.
with only one interior breakpoint. Starting with normalized accumulated chord length as the initial parameterization r(O), and the spline given in Fig. 7 (a) as s(O), we generate iterates r(l) and s(l), I = 1,. . ., from (5) and (6) as described in Section 2.1. Results for 1 = 10, 100, and 1000 are given in Fig. 8. The curve obtained after 10 iterations, Fig. 8 (a), already appears to be a better fit to the data than the cubic spline with three interior knots shown in Fig. 7 (c). In Fig. 8 (b), at I = 100, some additional improvement is evident with a noticably better lit near the peak. Continuing the iteration produces significant changes only near this peak. We were somewhat surprised by the pronounced overshoot that developes after 1000 iterations (Fig. 8 (c) ). This behavior suggests that the minimum norm solution may not be desired in all cases. Indeed, the application highlighted here is probably better served by stopping after 100 iterations. We also point out that the results in Fig. 8 (c) at 1 = 1000 are not yet optimal. Continued iteration, with 1 as large as 5000,
251
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
0.75
-0.25
(a)
0.75 (b)
0 iterations
1.75/
\
10 iterations
2.0
1 .o
1 .o
0.0 0.75
4 I.2 (c)
I
1.75
I
+
-0.25
100 iterations
0.75 (d)
I
I
-
1.75
,
’
*
800 iterations
Fig. 9. Example to illustrate nonexistence.
produces further slight reduction in the objective as well as a larger overshoot and somewhat sharper peak. This slow convergence (or non convergence) was also noted in (Plass and Stone, 1983), but not explained. The illustration which we provide below in Example 3 and our results in Section 3.1 suggest that the problem is caused by a lack of existence that arises for certain data. It is also important to note that the results shown in Fig. 8 are not sensitive to the precise placement of the single interior breakpoint. The experiment was repeated with the interior knot placed at 0.3, 0.4, 0.6 and 0.7. There were no significant visual differences among these cases. Similar results are also obtained using more than one interior breakpoint.
Example 3. The third example in this section question that we will discuss in Section 3.
will illustrate
the existence
258
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
3 60.0-z 3 I E 40.0Q) t E 8 20.0.
0.0
1 0
200
400
600
800
Iteration Number Fig. 10. B-spline coefficient
magnitude
vs. iteration number.
We start with data (Xi, yi ), i = 1, . . . ,33 where Xi =
0,
xi = 1, Xi =
2,
(i- 1) ---R-’ yi = (22 - ‘), yi =
yi =
(i ‘023) 10 ’
i = l,...,
11,
i = 12,. . . ,22, i = 23,...,33.
Again, we parametrize the data initially by accumulated chord length and attempt to fit it with a parametric cubic spline. This time we use cubic polynomials to approximate the components. To start calculations we compute the ordinary least squares fit. This is shown, along with the data, in Fig. 9(a). After 10 iterations the results displayed in Fig. 9(b) are obtained. At this point, the least square error is 0.057, while the maximum magnitudes of the B-spline coefficients for the x and y components are, respectively, 2.17 and 8.90. Similar results are shown after 100 and 800 iterations in Fig. 9 (cd). The behavior of the iterates in this example suggests that, while the error in approximation is decreasing to zero, the magnitude of at least one of the B-spline coeffficients associated with the y-coordinate of the fitted curve is growing without bound. This is confirmed graphically in Fig. 10, where we show the maximum B-spline coefficient magnitude plotted as a function of iteration number. While this example is somewhat pathological, the behavior that it illustrates can occur in more reasonable situations. Example 2 above illustrates this possibility. We also note that the choice of data for this example is more complicated than it need be. The same point could be made by choosing two, rather than three, lines of data points. We will return to a similar, but simpler, example in Section 3 to analytically illustrate the possibility of nonexistence.
Example 4. In the last example we illustrate the approximation
power of ODR
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 1I (1994) 247-267
259
1.0
ODR-
0.0
-1.c
-1
.o 1.c)
0.0
1.0
1.0’
0DR-e
s$Z)
-c
0.0
LSO
-1.0.
-1.0 -
1.0
-1.o
-1
.o
1.0
1.0.
1.0 ’
ODR -
-l.OU
-l.O(_,
-1.o
0.0
Fig. 11. ODR spline and ordinary
1.0
-1.o
least squares (LSQ) spline fits to logarithmic various cubic spline spaces.
1.0
spiral data for
splines. Before turning to a specific example, we make some brief comments related to this issue. Let F : [0, l] --) Rd be a smooth parametric curve in Rd. Additionally, let h = F(zi), i = l,...,N, where 0 < 71 <72<.‘.<7N< 1. Then the parametric spline s(q), with each component in Sf, determined as the ordinary least squares fit to the data (pi, f;:), i = 1,. . . , N satisfies the
S.P. Marin. P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
260
-16.0 I_, -5.0 -4.0
-3.0
-2.0
-1.0
0.0
In 00 Fig. 12. Example 4. Approximation
standard
estimate
error vs. knot spacing.
(see (de Boor, 1978 ) )
(10) where K spacing. good an the ODR
is a constant depending only on k and d, and ItI is maximum knot By its definition, the ODR spline, when it exists, will be at least as approximation. That is, if 7;, i = 1,. . . , N and s*, Pi (s*) E St, solve problem (2), then
(&s*~7;)lly2
< (~llJ-s(Ti),ly2 i=l
i=l 6
K&v
IIF(k)lloo
Itlk.
(11)
Our second point is the fact that, in general, the ODR spline will not improve upon the order of approximation. The effect of this approach is not to change the character of the approximation space, which might alter the order of approximation, but rather, to change the character of the data being approximated. Thus, the improvement comes as a result of reducing the contribution of IIF(k)llm to the error bound, with F the underlying interpolant of the data. As a specific illustration, we consider the parametric curve defined by
J’(v) =
emVcos (37~~) ePV sin (3ny)
>’
rc [O,ll.
As ‘1 ranges from 0 to 1, F (q ) traces out one and a half circuits of a logarithmic spiral. We obtain data by sampling the logarithmic spiral at 300 equally spaced points in [0, 1 ] . With this data we solve both the ordinary least
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
squares problem (3) and the ODR spline problem spline approximation spaces,
261
(2) on a sequence of cubic
Here, t(j) is a knot set consisting of four coincident knots at each of the endpoints -0.01 and 1.Ol, together with an additional 2-j - 1 uniformly spaced interior knots. Hence, at each refinement, we are halving the knot spacing. Representative geometric results are plotted in Fig. 11. On the left we have plotted, from top to bottom, the parametric cubic splines with components in $119 $2,) and S$,, computed using ordinary least squares. The plots on the right half show analogous results for the ODR splines. In each plot the solid curve is the fitted parametric spline, while the dashed curve is the logarithmic spiral given above. For smaller knot spacings, no discernable differences among the curves is present. For a more quantitative picture of the approximation characteristics in this example, we refer to Fig. 12. Here we plot In (a) as a function of In(h) where E is error and h is the knot spacing. Two curves are shown. The solid curve is the piecewise linear interpolant of the error data for the ordinary least squares results. The dashed curve is the piecewise linear interpolant of the error data for the ODR splines. As expected, the ODR results improve upon those obtained via ordinary least squares, but exhibit the same order of approximation for small h. For the coarser meshes, however, the improvement obtained can be significant. This observation points to perhaps the main motivation for considering ODR splines: To provide accurate, economical spline approximations to given geometric shapes. 2.3. Relationship to orthogonal distance regression The ODR spline problem considered here and in (Hoschek, 1988) is related, by way of their common objective of minimizing orthogonal distance, to the nonlinear orthogonal distance regression problem treated earlier in (Boggs et al., 1987), where the goal is to determine parameters /3 so that a function of the form Y = _f-(X,B)
minimizes the sum of squares of orthogonal distances from given data (xi, yi), i = I,... , n to the curve (or surface) defined by y = f (x, /3). Here, although x may be a vector, the dependent variable y is assumed to be a scalar, yielding an explicit representation for the fitted curve or surface. In this explicit ODR problem, the underlying optimization problem is unconstrained, and is solved in (Boggs et al., 1987) by a trust region method. Details of the problem’s structure are effectively used to formulate an efficient numerical solution procedure. The parametric problem we treat here is, in one respect, simpler. The component functionals are restricted to a spline space so that the unknown B-spline coefficients enter linearly in the definitions of these functionals. However, the parametric formulation of the problem introduces a vector dependent variable
262
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
and requires that inequality constraints be imposed on the components of the unknown vector r. Both the inequality constraints and the presence of a vector, rather than scalar, dependent variable change the character of the problem somewhat. Although, after accounting for these changes in structure, it may be possible to apply the methods of (Boggs et al., 1987), we have opted for the simpler implementation described in Section 2.1 in order to illustrate the utility of ODR splines. As noted in Section 1, the extension of the methods in (Boggs et al., 1987) to a parametric setting, as described in (Boggs et al., 1989), yields a somewhat different problem. The multiple response problem discussed in (Boggs et al., 1989) is again an unconstrained optimization problem. But, more significantly, its objective function includes a term that measures the difference between a candidate parameter value, ri, and a specified initial value, ri”‘. The implication of this is that, at optimality, errors are not orthogonal to the fitted curve at the determined parameter locations. The necessary conditions for an extrema in the setting of (Boggs et al., 1989) imply that the error vector makes an angle with the tangent to the fitted curve that depends on the initial parameter value r!‘). In the present approach, errors are orthogonal to the fitted curve at the optimized parameter values, as noted below in Section 3.4. This distinction might, in certain circumstances, lead to significantly different results.
3. Theoretical results This section further addresses the data fitting problem described in Section 1 and illustrated in Section 2 by considering questions of existence and uniqueness, and by providing a partial characterization of the solution to such problems. There are many open questions. One interesting general question which we will not touch on here is the following: What is the minimal degree parametric polynomial which interpolates the given data f? 3. I. Example of nonexistence In this subsection we present a simple example to show that under certain circumstances there will not exist a solution to the general ODR problem. Consider the situation where the parametric functions are quadratic splines with no interior knots (i.e. quadratic polynomials). We approximate the four data points fi = (O,O),
A = (O,l),
f3 = (l,l),
We will now show that the infimum s(7#:Pj(s)
“Ll = (l,O).
for the ODR problem
ES,3,0
<
71
<
‘*.
d
74
is zero. That is, <
= 0
1 >
S.P. Marin. P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
263
where t has only 6 knots (3 stacked at each end). This fact follows by observing that one can choose Pi(S)(X)
=x,
P2(s)(x)
=
(12)
4Px(l
(13)
-x),
where p is a parameter to be chosen later. This choice interpolates the data at 0 and 1. If we evaluate this parametric curve at 72 = 1/ (4p) and 73 = 1 - 1/(4p), we see that the error at these points goes to zero as j3 --+ 00. The computation is as follows:
IlShbh112 = =-.
($)2+(’-4($)(1 - +))’ 1
(14)
(15)
8P2
A similar computation will yield the same error at 73 since there is a symmetry in the data and function values about x = l/2. It follows that the error goes to zero as p goes to infinity. The above calculation shows that the ODR error for this problem is zero. It is easy to see that this error cannot be obtained by any selection of the 8 parameters (3 coefficients for each polynomial and the interior nodes 72 and 73). Thus we must have nonexistence in this simple (but not trivial) case. It is not clear in general which data sets will allow existence and which will not. This appears to be a very interesting and difficult problem. 3.2. Example of nonuniqueness In this subsection we present a simple nontrivial example to illustrate the possibility of multiple solutions to the ODR problem. The example we will present is due to Meyer (1990) and is a specialization of more general comments about multiplicity of solutions to linear least squares problems which Meyer relates in (Meyer, 1978 ) . The data which we will approximate are the three points fi = (-X&-c),
f2 = (0,2c),
f3 = (X&,-c)
where c is positive but otherwise arbitrary. These data are the vertices of an equilateral triangle whose centroid is at the origin (0,O). We will compute the ODR spline fit to this data assuming that the underlying spline space is the space of linear polynomials. Thus, the candidate parametric splines can be written in the form s(V) = (0,b)
+ ~(cos0,sin8),
--oo < q < 00,
(16)
where 19 is the angle the line makes with the positive x-axis and b is the y-intercept. To avoid ambiguity, we set b = 0 whenever 8 = n/2.
264
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
We find by direct calculation that the sum of squares distances measured from 3 to the line ( 16) is given by km?,,5
-s(v)]12
= 6c2 + 3b2cos28.
of the orthogonal
(17)
i=l
min,
Thus, Clcl llfi -s(r)ll 2 is minimized and let rT, i = 1,2,3 be the solutions of mF]]J;-s(r)]j2
i = 1,2,3
= I]~;:-.s(~T)]]~,
then r; < r; < r; holds whenever any line S(q) = q(cos6,sin8),
by choosing b = 0. If we fix b = 0
--~/6 < 8 < ~16. Thus, with 8 so restricted,
q E [ri,r;]
solves the ODR problem for this data. In (Meyer, 1978), Meyer has characterized those discrete point sets in [w2 which admit multiple linear ODR approximates. Interesting related, and to our knowledge, open questions remain regarding such issues for higher dimensional data and for higher degree polynomial approximates. 3.3. Remarks While the difficulties discussed in the previous two subsections might be seen as pessimistic, they can be easily handled numerically. One can recover existence by merely restricting the coefficients to satisfy some bound (or a bound on their differences). Nonuniqueness is generally not a problem since the minimization techniques which we recommend always lead to a smaller objective value. It is usually the case that a small objective value is the important criterion (as opposed to finding the smallest objective value).
3.4. Partial characterization To gain more insight into the parametric ODR spline problem described in (2), we derive the necessary conditions that characterize its solution. As noted in Section 2.1, we will study the form of the problem obtained by augmenting (2) with the order constraints 0 < ri < . . . < zm < 1. With this in mind, we consider the Lagrangian:
L(AC) :=
2 IlA i=l
-S(C,~i)l12
+
gii(ri-zi+l) i=O
where
S(C,X) :=
and ro := kCjBj,k,t(X) j=l
0, z,+~ := 1.
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
265
If there is a solution to the ODR problem (i.e., an optimal c* and r* with 0 < r: < r;,, G 1 for i = I,..+, m - 1) then the standard Kuhn-Tucker conditions apply yielding a A > 0 satisfying 0 = -2(f;-s(c*,r~),(s*)‘(r~))
i = l,...,m
+&-A&i,
i =
0 = &(TT - TT,,),
. . , m.
Here, (s*)‘(x) denotes the derivative of s(c*,x) with respect to x. Three interesting results immediately arise from studying the Lagrangian. First, if r;_i < 7T < 7r+1, then [J - s*(c*,t;)]
(s*)‘(z:),
I
which is the usual orthogonal regression condition. Second, if p of the nodes in r* coincide in the interior of [ 0, 1 ] (z:_ 1 < z; = . . . = T;+~_~ < T;+~), then
CTiAA+j Third, 7; <
77,
s(c*,
P
[
(s*)‘(zf).
I 1
if we have p > 1 nodes at the left endpoint then
51 = 0 (0 = r; = . . . =
7;+1),
0.-g=
cp=, “6
(
P
-
s(c*,z;),
(s*)‘(z;)
If we have p nodes at the right endpoint 7; = I), then
o+
CTihh-j
(
P
.
1
-
s(c*,zt,),
zm
=
1
(z&_~
(s’)‘($)
<
7&p+l
=
. . . =
.
)
There is another point to make. If the spline space we are working with reduces to a polynomial subspace (due to the fact that we have no interior knots), then it is easy to see that we must have orthogonality at the end points. That is, if p nodes coincide with 71 = 0, then
ZTIhJ;:+j-s(c*,zT) [
P
1 I
(s*)‘(sT),
i = 1.
Since polynomial subspaces are invariant under a linear change of variables, we see that the constraints 71 >, 0 and 7 ,,, < 1 are not binding, hence the result. There is a similar statement for the right end of the interval.
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
266
4. Approximation of functions While we have been motivated by the approximation of discrete data, we must comment on the connection ODR splines have to approximation of functions. Suppose we are given a curve F in Rd. Then the classical measure of discrepancy between F and an approximation s is 1
s
IIF
- s(dl12dT.
(18)
0
This measure may be inappropriate in certain situations. For example, if only the shape of the curve F is important (that is, if the speed of traversal of the curve is unimportant) then one should consider the error measure
E(s, F) := itf { ],,F(T)
-s(n(i)),l’d~},
(19)
0
where (Yis an increasing homeomorphism We briefly consider the problem
of [0, 1 ] into itself.
inf{E (s, F) : Pi (s) E SF}
(20)
Under appropriate hypotheses of smoothness one can derive upper bounds on the error from the classical spline error bounds as follows. First notice that if Q is a diffeomorphism, we can make the change of variables 7 = Q (z) yielding
Thus if l/a’ < COwe can apply classical error estimates inf{(E(s,J’))
‘I2 : Pi(s) E S,“} < Kinf $ ,II(F a II I/
where K is a constant
depending
yielding 0 Q-‘)(~)JJ~ \tlk,
(22)
only on k.
5. Summary We hope that the ideas presented here provoke more discussion and experimentation. The general problem of recovering a geometric curve (or surface) when no explicit parameterization is given appears to be a pervasive problem. While we are unable to prove that ODR splines have favorable shape preserving properties, our numerical experiments suggest that they can play a very useful role in faithfully representing geometric shapes with fewer parameters than are generally required in competing methods.
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
261
Acknowledgement We would like to thank Weston Meyer for his insightful comments regarding uniqueness for ODR splines, and for suggesting an alternate, more direct counterexample to uniqueness.
References Boggs, P.T., Byrd, R.H. and Schnabel, R.B. (1987), A stable and efficient algorithm for nonlinear orthogonal distance regression, SIAM J. Sci. Stat. Comput. 8, 1052-1078. Boggs, P.T., Donaldson, J.T., Byrd, R.H. and Schnabel, R.B. (1989), ODRPACK: Software for weighted orthogonal distance regression, ACM Trans. Math. Software 15 (4), 348-364. de Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Foley, T.A. and Nielson, GM. (1989), Knot selection for parametric spline interpolation, in Lyche, T. and Schumaker, L.L., eds., Mathematical Methods in Computer Aided Geometric Design, Academic Press, Boston, MA, 261-267. Herz, R.K. and Marin, S.P., (1980), Surface chemistry models of carbon monoxide oxidation on supported platinum catalysts, J. Catalysis 65, 281-296. Hoschek, J. (1988), Intrinsic parameterization for approximation, Computer Aided Geometric Design 5, 27-31. Lee, E.Y.T. ( 1989), Choosing nodes in parametric curve interpolation, Computer-Aided Design 21, 363-370. Marin, S.P. (1984), An approach to data parametrization in parametric cubic spline interpolation problems, J. Approx. Theory 41, 64-86. Meyer, W.W. (1978), Orthogonal projections. Solution to problem 988. Math. Magazine 51, 71-72. Meyer, W.W. (1990), Private communication. O’Dell, C.L. (1985), Approximating data with parametric B-splines for computer aided design, M.S. Thesis, Department of Computer Science, University of Utah, Salt Lake City. Pinkus, A. (1988), On the smoothest interpolation, SIAM J. Math. Anal. 19, 1431-1441. Plass, M. and Stone, M. (1983), Curve-fitting with piecewise parametric cubits, Computer Graphics 17, 229-239. Rademacher, C. and Scherer, K. (1990), Algorithms for computing best parametric cubic interpolation, in: Mason, J.C. and Cox, M.G., eds., Algorithms for Approximation II, Chapman and Hall, London, 193-207. Scherer, K. (1987), Uniqueness of best parametric interpolation by cubic spline curves, preprint. Scherer, K. and Smith, P.W. (1989), Existence of best parametric interpolation by curves, SIAM J. Math. Anal. 20, 160-168. Toepfer, H.J. (1981), Models for smooth curve fitting, Numererische Methoden der Approximationstheorie, 6, Internat. Ser. Numer. Math., Vol. 59, Birkhauser, Basel.
Computer Aided Geometric Design 11 (1994) 247-267
Parametric
approximation Samuel
of data using ODR splines
P. Marin a,*, Philip
W. Smith b
a Mathematics Department, GM Research and Development Center, Warren, MI 48090-9055, USA b 14311 Kellywood, Houston, TX 77079, USA
Received June 1990; revised August 1991
Abstract We consider curve (that is, and numerical arises naturally
the problem of approximating data in space that arises from a geometric a curve with no explicit parameterization given). We treat theoretical aspects of a parametric orthogonal distance regression problem which in the context of geometric data fitting.
Key words: Orthogonal B-spline curves
distance
regression;
Parametric
approximation;
Least squares;
1. Introduction Given m ordered data points in Rd, {A}, i = 1,. . . , m, we want to find a curve which approximates the data and has a relatively simple form. By a curve we mean a function g which maps the unit interval continuously into Rd. Consider, for example, the data shown in Fig. 1. The points here were measured using a coordinate measuring machine along a cross section of a spiral surface of an air compressor component. The goal is to fit a parametric curve to the data and thus analytically represent the shape of the cross section. We will measure the approximation power of curves fitted to such data by computing E(g,f):=min 1
eJ/f;-g(Ti)J/‘:~iE it
[O,l],i=
l,...,m}.
(1)
I
We emphasize that this is not the classical measure of discrepancy which would require us to choose a priori the points where the errors are to be measured. Of course, if no restriction is placed on g, then we could just interpolate the *Corresponding author. Email: [email protected]. 0167-8396/94/$07.00 @ 1994 Elsevier Science B.V. All rights reserved 0167-8396(93)E0029-D
SSDI
248
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
D
Y (mm) 0.0 0
D
-50.0 _j
1
-50.0
0.0
50.0
x (mm) Fig. 1. Sample CMM data for parametric
curve in R2.
data at ri = (i - 1 )/ (m - 1). This would be unacceptable in most situations where the data are noisy. Furthermore, it would mean that the function g would become more complex as m grows. We avoid this issue of complexity by restricting the coordinate functionals of g to lie in a finite dimensional spline subspace. This leads us to the following fundamental approximation problem which we consider in this paper: inf{E(s,f):Pi(s)ES~,i=
l,...,
d}.
(2)
Here, SF is the subspace of kth order spline functions with knot sequence t (see de Boor (1978)) and Pi is the projection operator onto the ith coordinate. Thus s is a curve in Rd with each coordinate functional a spline in Sf. In the following we assume that both k and t are known and fixed. We also require that the interval spanned by the knot sequence t contains [0,1 ] so that the spline functions are defined over the target parameterization interval. We call solutions to (2) (if they exist) ODR splines. The ODR stands for orthogonal distance regression. This is discussed in (Boggs et al., 1987) where an algorithm is presented and analyzed to solve the ODR problem for explicitly represented curves and surfaces. The parametric problem posed and studied in our paper represents a variant which is not directly addressed by the methods outlined in (Boggs et al., 1987). In (Boggs et al., 1989), an extension of the methods in (Boggs et al., 1987) to multiple response data is described. This provides a capability to handle parametric problems. However, the objective function used in (Boggs et al., 1989 ) differs from (2) and therefore leads to a different geometric interpretation of the results. We say more about
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
-50.0 _! -50.0
249
0.0 x (mm)
Fig. 2. Ordinary least squares fit to data from Fig. 1.
this in Section 2.3. The parametric problem we consider has been previously posed by Plass and Stone ( 1983) and by Hoschek ( 1988) and references given there. These references present numerical approaches together with some examples from CAD applications. Additional examples and comparisons to other approaches are given by O’Dell ( 1985 ). We provide additional insights pertaining to the existence, uniqueness, characterization, and the approximation power of the solutions to these problems. A core element in the parametric ODR problem is that of data parameterization. For related literature covering other aspects of data parameterization for parametric spline curves we refer to (Foley and Nielson, 1989; Lee, 1989; Marin, 1984; Pinkus, 1988; Rademacher and Scherer, 1990; Scherer and Smith, 1989; Scherer, 1987; Toepfer, 198 1). There is a clear difference between the approximation problem above and the usual method for solving the parametric curve fitting problem. The standard approach to the data fitting problem requires the selection of a subspace as above in which the coordinate functionals must lie. Additionally, points 0 6 rr d *.. d rm Q 1 are selected. Now the standard least squares problem is solved (3) This standard least squares approach may be applied to the data of Fig. 1 by first choosing, for example, {7i, i = 1, . . . , m} to be normalized accumulated chord length and SF to be piecewise cubits defined on [0, 1 ] with seven
250
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
Y (mm)
0.0 -
-50.0 J -50.0
0.0
50.0
x (mm) Fig. 3. ODR spline fit to data from Fig. 1.
uniformly spaced interior breakpoints. The results are shown in Fig. 2. This standard least squares problem is linear. The price of linearity, however, is the inflexibility of prescribing the points ri. One of the principal advantages of the ODR spline problem (2) is that it incorporates the selection of {r} as part of its solution and thus alleviates this inflexibility. To see what might be gained by resorting to an ODR spline fit, we present, in Fig. 3, the parametric, piecewise cubic solution of (2) for the data shown in Fig. 1. The spline space SF is the same here as in the ordinary least squares case discussed above. These results, however, show a noticeable improvement over those displayed in Fig. 2. In the remainder of the paper we expand our discussion of this approach to parametric curve fitting according to the following outline. In Section 2 we describe a solution algorithm and present several numerical examples which provide additional illustrations of the capabilities and limitations of ODR splines. Section 3 contains theoretical results addressing questions of existence, uniqueness, and the characterization of solutions of these problems. In Section 4 we comment on underlying approximation issues. We conclude with a brief summary in Section 5.
2. Numerical examples The main purpose of this section is to present several examples which illustrate both the pitfalls and the power of ODR splines versus the classical
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
251
least squares method. Solutions for the examples chosen are obtained using a simple iterative technique. The method, which for completeness we describe below, is essentially that presented in (Hoschek, 1988). At the end of this section we comment on the connection between the data fitting problem as we have described it and orthogonal distance regression. 2.1. Algorithm description The starting point is the observation that if s* (q) solves (2), then it is necessary that s* also solve the ordinary least squares problem (3) with ri, i = l,..., m chosen according to
s*(7i)ll
Il_h-
= ,@, 1l.h- s*(rl)lI. . .
(4)
The algorithm described in (Hoschek, 1988) proceeds by iterating between problems (3) and (4). Starting with an initial parameterization {r(O)} and a spline s(O), (s(O) solves (3) with ri = ri”, i = 1,. . .,m), s(I) are generated which satisfy, respectively,
iterates
{rj”}
and
and
min
{
2 IIf;: - scr)(711,)l12 : Pj(s(l))
E
s:).
(6)
i=l
In the implementation discussed in (Hoschek, 1988)) the univariate minimization problem (5 ) is solved approximately by first constructing, at q = $-I) a local linear approximation to the parametric curve s (I-l ) (v ). With ,11-l,‘. m (5) replaced by its linear approximate, an updated parameter value, ,!I) from I ) is determined
7% I
(S+lqp),fi =
p
+
v
[
- su-l)(p))
~I,cr-w(7!~-9~~2 I
1.
Here, v is a scalar search variable, initially equal to 1, which is halved until the predicted rj” satisfies [IA - s”-‘)(s~“)ll < llfi - ~(~-‘)(r~‘-“)ll. In most cases we used this approach. However, in a few demanding examples where the fitted curve developed regions of high curvature, this local method was not appropriate. To proceed in these situations, we relied on a simple global, uniform search to solve (5 ) . Normally, we expect that the fitted curve will preserve the order of the data implied by the initial parametrization. That is, if r;, i = 1,. . . , m is the final parametrization determined by solving (2)) then we desire that
252
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
It is clear that the minimization problem (5) may not yield solutions r!l’, i = 1, . . . , m which satisfy the monotonicity condition (7). The likeliness of this happening depends, in part, on the choice of the spline subspace. For example, if we chose to solve the ODR problem for the data of Fig. 1 using parametric splines whose projections are linear polynomials, then solutions of (5) would not form a monotone sequence as desired. The remedy here is to increase the order k and refine the knot sequence t so that the spline subspace 5’: can more closely approximate the data. Our limited experience suggests that, for large enough k and fine enough t, the monotonicity condition (7) will be satistied as the iteration proceeds. We can, however, envision pathological examples where this will not be the case. In some cases it may be desirable to include the constraints (7) explicitly in the problem formulation. It may also be possible to solve such order constrained ODR spline problems using an iterative technique similar to that described above for the problem without order constraints. To do this we would change the update procedure for the vector r (0 of parameter values. Instead of solving the sequence of univariate minimization problems given in (5 ) to obtain r!” ,I = l,..., m, these values would be derived as the solution of the multivariate problem
The iteration would proceed by cycling between (8 ) and (6 ) . Although we have not carried out numerical experiments to test this approach, we do derive, in Section 3.4, the necessary conditions for a solution in this setting. Without a guarantee of existence of solutions for this parametric ODR problem (see Section 3), it is difficult to comment extensively on the convergence of the numerical algorithms. One useful point can be made, however. The iterates {r(“} and {s(l)} constructed according to (5 ) and (6 ) can be shown to satisfy E(s(‘-‘),f)
2 E(s’“,f),
1 = l)...
(9)
E (s, f ) is defined by ( 1). Thus the sequence of objective values {E (s(l), f )} is monotone decreasing, bounded below and hence converges to a value E*. As we will illustrate later, the spline coefficients associated with
where
the iterates {.s(‘)} may, however, grow without bound. We now consider several examples to highlight the capabilities limitations of ODR splines.
and the
2.2. Examples In the examples that follow, we seek ODR spline fits to four two-dimensional problems. In each case, the data is initially parameterized on the unit interval [0,1 ] and the underlying spline space $ is chosen so that the interval spanned by its knots strictly contains [0, 11. The only constraints imposed in solving the
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
0.5
I/
0.5
0.0
1.0
253
0.5
0.0
1.0
7
7
Fig. 4. Ordinary least squares fit to circle data.
Fig. 5. ODR spline fit to circle data.
minimization problem (5) are that 0 d ~j” d 1, i = 1,. . . , m. The graphs that accompany these examples show the parametric ODR spline, s*(q), plotted for q E [r;,z&].
Example 1. For our first example we construct the solution of problem (2) when the data {J;:}, i = 1,. . . , m, consist of 30 points in [w2uniformly spaced on the unit circle. The spline space consists of cubic B-splines with 11 knots t1 = t2 = t3 = t4 = -0.01, t5 = 0.245,
t(j = 0.5,
t8 = t9 = tlo = tll
= 1.01.
t7 = 0.755,
254
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
-ll.Ol
.
-10.0
. 1 -9.0
,
,
, , -8.6
,
,
( -7.0
109 (Cco)
Fig. 6. Piecewise linear interpolant of sample reaction rate data.
The data is parameterized (0)
T1
=
0,
,y
=
initially by normalized
T!O)+ r-1
IN - “LIII Ci”=,llfi
-jj_lI1’
accumulated
chord length,
.
E= 2,***,m*
The ordinary least squares fit to this data is given in Fig. 4. A plot of the fitted parametric spline is shown along with a plot of curvature. The same information is given in Fig. 5, but for the ODR spline fit to the data. Notable differences between Figs. 4 and 5 concern the behavior of the parametric curves near the data point ( 1,O) (no periodicity condition was imposed in either case) and the magnitude of the observed oscillations in the curvature plots. These are significantly lower in the case of the ODR spline. We note that periodicity can be imposed in the spline subspace without changing the solution methodology. We chose not to do so in this example to illustrate how well the base method captures geometric properties of the data.
Example 2. The second example
illustrates a potential application of these methods to laboratory data fitting. The situation we have in mind concerns a reacting chemical system which may possess multiple steady states. To obtain data for the example, we have sampled one of the model curves given in (Herz and Marin, 1980). The abscissa, Cco, of each data point is an input concentration of carbon monoxide in a catalytic system. The ordinate, R, is the steady state oxidation rate achieved by the system. In Fig. 6 we have plotted (log (Cc0 ),log (R ) ) at each sampled point and, to establish order, the piecewise linear interpolant of the data. Since multiple reaction rates are attained at a single input concentration, parametric curve fitting methods must be used.
255
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
-8.0
0
-8.0
OS0
(b)
(a) -9.c
-9.0
lo9 w
lo9 (W
-10.0
!
-11.o
-10.0
-8.0
-9.0
-7.0
-10.0
-9.0
-8.0
-7.0
log (Cco)
b2(Cco)
-8.0
-9.0
log (R) -10.0
-11.0 I_. -10.0
-9.0
-8.0
-7.0
~og~cco)
Fig. 7. Ordinary
least squares fits to reaction rate data. (a) One interior knot. (b) Two interior knots. (c) Three interior knots.
With the data parameterized by normalized accumulated chord length we show, in Fig. 7 (a), the ordinary least squares fit using parametric cubic splines having breakpoints at -0.01, 0.50, and 1.Ol. With only one interior breakpoint, the data are not well fit. One way to improve the fit is to introduce more breakpoints. In Figs. 7(b) and 7(c) we show the ordinary least squares fit obtained using, respectively, two and three uniformly spaced interior breakpoints in the interval (-0.01, 1.01). The curve in Fig. 7 (b) is still a rather poor fit to the data. It can be improved somewhat by more carefully positioning the two interior breakpoints. The curve shown in Fig. 7 (c) (three interior knots) is better, but still fails to model the behavior near the peak in the data. The ODR spline fitting technique offers an alternative procedure that is not so sensitive to the number or location of breakpoints. In fact, we can very effectively fit the reaction rate data with a parametric cubic spline curve
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
256
-8.0
-8.0
1
-9.0
-9.0.
‘09 CR)
‘og CR)
-10.0.
-10.0.
-1l.OI
-11.0I_
-10.0
1
-9.0
-8.0
-7.0
-10.0
109 (Cc,)
-9.0
-8.0
-7.0
~og(cco)
-8.0
-11.0A -10.0
-9.0
-8.0
-7.0
'09Gzo)
Fig. 8. ODR fit to reaction
rate data (one interior knot). (c) 1000 iterations.
(a) 10 iterations.
(b) 100 iterations.
with only one interior breakpoint. Starting with normalized accumulated chord length as the initial parameterization r(O), and the spline given in Fig. 7 (a) as s(O), we generate iterates r(l) and s(l), I = 1,. . ., from (5) and (6) as described in Section 2.1. Results for 1 = 10, 100, and 1000 are given in Fig. 8. The curve obtained after 10 iterations, Fig. 8 (a), already appears to be a better fit to the data than the cubic spline with three interior knots shown in Fig. 7 (c). In Fig. 8 (b), at I = 100, some additional improvement is evident with a noticably better lit near the peak. Continuing the iteration produces significant changes only near this peak. We were somewhat surprised by the pronounced overshoot that developes after 1000 iterations (Fig. 8 (c) ). This behavior suggests that the minimum norm solution may not be desired in all cases. Indeed, the application highlighted here is probably better served by stopping after 100 iterations. We also point out that the results in Fig. 8 (c) at 1 = 1000 are not yet optimal. Continued iteration, with 1 as large as 5000,
251
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
0.75
-0.25
(a)
0.75 (b)
0 iterations
1.75/
\
10 iterations
2.0
1 .o
1 .o
0.0 0.75
4 I.2 (c)
I
1.75
I
+
-0.25
100 iterations
0.75 (d)
I
I
-
1.75
,
’
*
800 iterations
Fig. 9. Example to illustrate nonexistence.
produces further slight reduction in the objective as well as a larger overshoot and somewhat sharper peak. This slow convergence (or non convergence) was also noted in (Plass and Stone, 1983), but not explained. The illustration which we provide below in Example 3 and our results in Section 3.1 suggest that the problem is caused by a lack of existence that arises for certain data. It is also important to note that the results shown in Fig. 8 are not sensitive to the precise placement of the single interior breakpoint. The experiment was repeated with the interior knot placed at 0.3, 0.4, 0.6 and 0.7. There were no significant visual differences among these cases. Similar results are also obtained using more than one interior breakpoint.
Example 3. The third example in this section question that we will discuss in Section 3.
will illustrate
the existence
258
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
3 60.0-z 3 I E 40.0Q) t E 8 20.0.
0.0
1 0
200
400
600
800
Iteration Number Fig. 10. B-spline coefficient
magnitude
vs. iteration number.
We start with data (Xi, yi ), i = 1, . . . ,33 where Xi =
0,
xi = 1, Xi =
2,
(i- 1) ---R-’ yi = (22 - ‘), yi =
yi =
(i ‘023) 10 ’
i = l,...,
11,
i = 12,. . . ,22, i = 23,...,33.
Again, we parametrize the data initially by accumulated chord length and attempt to fit it with a parametric cubic spline. This time we use cubic polynomials to approximate the components. To start calculations we compute the ordinary least squares fit. This is shown, along with the data, in Fig. 9(a). After 10 iterations the results displayed in Fig. 9(b) are obtained. At this point, the least square error is 0.057, while the maximum magnitudes of the B-spline coefficients for the x and y components are, respectively, 2.17 and 8.90. Similar results are shown after 100 and 800 iterations in Fig. 9 (cd). The behavior of the iterates in this example suggests that, while the error in approximation is decreasing to zero, the magnitude of at least one of the B-spline coeffficients associated with the y-coordinate of the fitted curve is growing without bound. This is confirmed graphically in Fig. 10, where we show the maximum B-spline coefficient magnitude plotted as a function of iteration number. While this example is somewhat pathological, the behavior that it illustrates can occur in more reasonable situations. Example 2 above illustrates this possibility. We also note that the choice of data for this example is more complicated than it need be. The same point could be made by choosing two, rather than three, lines of data points. We will return to a similar, but simpler, example in Section 3 to analytically illustrate the possibility of nonexistence.
Example 4. In the last example we illustrate the approximation
power of ODR
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 1I (1994) 247-267
259
1.0
ODR-
0.0
-1.c
-1
.o 1.c)
0.0
1.0
1.0’
0DR-e
s$Z)
-c
0.0
LSO
-1.0.
-1.0 -
1.0
-1.o
-1
.o
1.0
1.0.
1.0 ’
ODR -
-l.OU
-l.O(_,
-1.o
0.0
Fig. 11. ODR spline and ordinary
1.0
-1.o
least squares (LSQ) spline fits to logarithmic various cubic spline spaces.
1.0
spiral data for
splines. Before turning to a specific example, we make some brief comments related to this issue. Let F : [0, l] --) Rd be a smooth parametric curve in Rd. Additionally, let h = F(zi), i = l,...,N, where 0 < 71 <72<.‘.<7N< 1. Then the parametric spline s(q), with each component in Sf, determined as the ordinary least squares fit to the data (pi, f;:), i = 1,. . . , N satisfies the
S.P. Marin. P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
260
-16.0 I_, -5.0 -4.0
-3.0
-2.0
-1.0
0.0
In 00 Fig. 12. Example 4. Approximation
standard
estimate
error vs. knot spacing.
(see (de Boor, 1978 ) )
(10) where K spacing. good an the ODR
is a constant depending only on k and d, and ItI is maximum knot By its definition, the ODR spline, when it exists, will be at least as approximation. That is, if 7;, i = 1,. . . , N and s*, Pi (s*) E St, solve problem (2), then
(&s*~7;)lly2
< (~llJ-s(Ti),ly2 i=l
i=l 6
K&v
IIF(k)lloo
Itlk.
(11)
Our second point is the fact that, in general, the ODR spline will not improve upon the order of approximation. The effect of this approach is not to change the character of the approximation space, which might alter the order of approximation, but rather, to change the character of the data being approximated. Thus, the improvement comes as a result of reducing the contribution of IIF(k)llm to the error bound, with F the underlying interpolant of the data. As a specific illustration, we consider the parametric curve defined by
J’(v) =
emVcos (37~~) ePV sin (3ny)
>’
rc [O,ll.
As ‘1 ranges from 0 to 1, F (q ) traces out one and a half circuits of a logarithmic spiral. We obtain data by sampling the logarithmic spiral at 300 equally spaced points in [0, 1 ] . With this data we solve both the ordinary least
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
squares problem (3) and the ODR spline problem spline approximation spaces,
261
(2) on a sequence of cubic
Here, t(j) is a knot set consisting of four coincident knots at each of the endpoints -0.01 and 1.Ol, together with an additional 2-j - 1 uniformly spaced interior knots. Hence, at each refinement, we are halving the knot spacing. Representative geometric results are plotted in Fig. 11. On the left we have plotted, from top to bottom, the parametric cubic splines with components in $119 $2,) and S$,, computed using ordinary least squares. The plots on the right half show analogous results for the ODR splines. In each plot the solid curve is the fitted parametric spline, while the dashed curve is the logarithmic spiral given above. For smaller knot spacings, no discernable differences among the curves is present. For a more quantitative picture of the approximation characteristics in this example, we refer to Fig. 12. Here we plot In (a) as a function of In(h) where E is error and h is the knot spacing. Two curves are shown. The solid curve is the piecewise linear interpolant of the error data for the ordinary least squares results. The dashed curve is the piecewise linear interpolant of the error data for the ODR splines. As expected, the ODR results improve upon those obtained via ordinary least squares, but exhibit the same order of approximation for small h. For the coarser meshes, however, the improvement obtained can be significant. This observation points to perhaps the main motivation for considering ODR splines: To provide accurate, economical spline approximations to given geometric shapes. 2.3. Relationship to orthogonal distance regression The ODR spline problem considered here and in (Hoschek, 1988) is related, by way of their common objective of minimizing orthogonal distance, to the nonlinear orthogonal distance regression problem treated earlier in (Boggs et al., 1987), where the goal is to determine parameters /3 so that a function of the form Y = _f-(X,B)
minimizes the sum of squares of orthogonal distances from given data (xi, yi), i = I,... , n to the curve (or surface) defined by y = f (x, /3). Here, although x may be a vector, the dependent variable y is assumed to be a scalar, yielding an explicit representation for the fitted curve or surface. In this explicit ODR problem, the underlying optimization problem is unconstrained, and is solved in (Boggs et al., 1987) by a trust region method. Details of the problem’s structure are effectively used to formulate an efficient numerical solution procedure. The parametric problem we treat here is, in one respect, simpler. The component functionals are restricted to a spline space so that the unknown B-spline coefficients enter linearly in the definitions of these functionals. However, the parametric formulation of the problem introduces a vector dependent variable
262
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
and requires that inequality constraints be imposed on the components of the unknown vector r. Both the inequality constraints and the presence of a vector, rather than scalar, dependent variable change the character of the problem somewhat. Although, after accounting for these changes in structure, it may be possible to apply the methods of (Boggs et al., 1987), we have opted for the simpler implementation described in Section 2.1 in order to illustrate the utility of ODR splines. As noted in Section 1, the extension of the methods in (Boggs et al., 1987) to a parametric setting, as described in (Boggs et al., 1989), yields a somewhat different problem. The multiple response problem discussed in (Boggs et al., 1989) is again an unconstrained optimization problem. But, more significantly, its objective function includes a term that measures the difference between a candidate parameter value, ri, and a specified initial value, ri”‘. The implication of this is that, at optimality, errors are not orthogonal to the fitted curve at the determined parameter locations. The necessary conditions for an extrema in the setting of (Boggs et al., 1989) imply that the error vector makes an angle with the tangent to the fitted curve that depends on the initial parameter value r!‘). In the present approach, errors are orthogonal to the fitted curve at the optimized parameter values, as noted below in Section 3.4. This distinction might, in certain circumstances, lead to significantly different results.
3. Theoretical results This section further addresses the data fitting problem described in Section 1 and illustrated in Section 2 by considering questions of existence and uniqueness, and by providing a partial characterization of the solution to such problems. There are many open questions. One interesting general question which we will not touch on here is the following: What is the minimal degree parametric polynomial which interpolates the given data f? 3. I. Example of nonexistence In this subsection we present a simple example to show that under certain circumstances there will not exist a solution to the general ODR problem. Consider the situation where the parametric functions are quadratic splines with no interior knots (i.e. quadratic polynomials). We approximate the four data points fi = (O,O),
A = (O,l),
f3 = (l,l),
We will now show that the infimum s(7#:Pj(s)
“Ll = (l,O).
for the ODR problem
ES,3,0
<
71
<
‘*.
d
74
is zero. That is, <
= 0
1 >
S.P. Marin. P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
263
where t has only 6 knots (3 stacked at each end). This fact follows by observing that one can choose Pi(S)(X)
=x,
P2(s)(x)
=
(12)
4Px(l
(13)
-x),
where p is a parameter to be chosen later. This choice interpolates the data at 0 and 1. If we evaluate this parametric curve at 72 = 1/ (4p) and 73 = 1 - 1/(4p), we see that the error at these points goes to zero as j3 --+ 00. The computation is as follows:
IlShbh112 = =-.
($)2+(’-4($)(1 - +))’ 1
(14)
(15)
8P2
A similar computation will yield the same error at 73 since there is a symmetry in the data and function values about x = l/2. It follows that the error goes to zero as p goes to infinity. The above calculation shows that the ODR error for this problem is zero. It is easy to see that this error cannot be obtained by any selection of the 8 parameters (3 coefficients for each polynomial and the interior nodes 72 and 73). Thus we must have nonexistence in this simple (but not trivial) case. It is not clear in general which data sets will allow existence and which will not. This appears to be a very interesting and difficult problem. 3.2. Example of nonuniqueness In this subsection we present a simple nontrivial example to illustrate the possibility of multiple solutions to the ODR problem. The example we will present is due to Meyer (1990) and is a specialization of more general comments about multiplicity of solutions to linear least squares problems which Meyer relates in (Meyer, 1978 ) . The data which we will approximate are the three points fi = (-X&-c),
f2 = (0,2c),
f3 = (X&,-c)
where c is positive but otherwise arbitrary. These data are the vertices of an equilateral triangle whose centroid is at the origin (0,O). We will compute the ODR spline fit to this data assuming that the underlying spline space is the space of linear polynomials. Thus, the candidate parametric splines can be written in the form s(V) = (0,b)
+ ~(cos0,sin8),
--oo < q < 00,
(16)
where 19 is the angle the line makes with the positive x-axis and b is the y-intercept. To avoid ambiguity, we set b = 0 whenever 8 = n/2.
264
S.P. Marin, P. W. Smith / Computer Aided Geometric Design I1 (1994) 247-267
We find by direct calculation that the sum of squares distances measured from 3 to the line ( 16) is given by km?,,5
-s(v)]12
= 6c2 + 3b2cos28.
of the orthogonal
(17)
i=l
min,
Thus, Clcl llfi -s(r)ll 2 is minimized and let rT, i = 1,2,3 be the solutions of mF]]J;-s(r)]j2
i = 1,2,3
= I]~;:-.s(~T)]]~,
then r; < r; < r; holds whenever any line S(q) = q(cos6,sin8),
by choosing b = 0. If we fix b = 0
--~/6 < 8 < ~16. Thus, with 8 so restricted,
q E [ri,r;]
solves the ODR problem for this data. In (Meyer, 1978), Meyer has characterized those discrete point sets in [w2 which admit multiple linear ODR approximates. Interesting related, and to our knowledge, open questions remain regarding such issues for higher dimensional data and for higher degree polynomial approximates. 3.3. Remarks While the difficulties discussed in the previous two subsections might be seen as pessimistic, they can be easily handled numerically. One can recover existence by merely restricting the coefficients to satisfy some bound (or a bound on their differences). Nonuniqueness is generally not a problem since the minimization techniques which we recommend always lead to a smaller objective value. It is usually the case that a small objective value is the important criterion (as opposed to finding the smallest objective value).
3.4. Partial characterization To gain more insight into the parametric ODR spline problem described in (2), we derive the necessary conditions that characterize its solution. As noted in Section 2.1, we will study the form of the problem obtained by augmenting (2) with the order constraints 0 < ri < . . . < zm < 1. With this in mind, we consider the Lagrangian:
L(AC) :=
2 IlA i=l
-S(C,~i)l12
+
gii(ri-zi+l) i=O
where
S(C,X) :=
and ro := kCjBj,k,t(X) j=l
0, z,+~ := 1.
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
265
If there is a solution to the ODR problem (i.e., an optimal c* and r* with 0 < r: < r;,, G 1 for i = I,..+, m - 1) then the standard Kuhn-Tucker conditions apply yielding a A > 0 satisfying 0 = -2(f;-s(c*,r~),(s*)‘(r~))
i = l,...,m
+&-A&i,
i =
0 = &(TT - TT,,),
. . , m.
Here, (s*)‘(x) denotes the derivative of s(c*,x) with respect to x. Three interesting results immediately arise from studying the Lagrangian. First, if r;_i < 7T < 7r+1, then [J - s*(c*,t;)]
(s*)‘(z:),
I
which is the usual orthogonal regression condition. Second, if p of the nodes in r* coincide in the interior of [ 0, 1 ] (z:_ 1 < z; = . . . = T;+~_~ < T;+~), then
CTiAA+j Third, 7; <
77,
s(c*,
P
[
(s*)‘(zf).
I 1
if we have p > 1 nodes at the left endpoint then
51 = 0 (0 = r; = . . . =
7;+1),
0.-g=
cp=, “6
(
P
-
s(c*,z;),
(s*)‘(z;)
If we have p nodes at the right endpoint 7; = I), then
o+
CTihh-j
(
P
.
1
-
s(c*,zt,),
zm
=
1
(z&_~
(s’)‘($)
<
7&p+l
=
. . . =
.
)
There is another point to make. If the spline space we are working with reduces to a polynomial subspace (due to the fact that we have no interior knots), then it is easy to see that we must have orthogonality at the end points. That is, if p nodes coincide with 71 = 0, then
ZTIhJ;:+j-s(c*,zT) [
P
1 I
(s*)‘(sT),
i = 1.
Since polynomial subspaces are invariant under a linear change of variables, we see that the constraints 71 >, 0 and 7 ,,, < 1 are not binding, hence the result. There is a similar statement for the right end of the interval.
S.P. Marin, P. W. Smith / Computer Aided Geometric Design II (1994) 247-267
266
4. Approximation of functions While we have been motivated by the approximation of discrete data, we must comment on the connection ODR splines have to approximation of functions. Suppose we are given a curve F in Rd. Then the classical measure of discrepancy between F and an approximation s is 1
s
IIF
- s(dl12dT.
(18)
0
This measure may be inappropriate in certain situations. For example, if only the shape of the curve F is important (that is, if the speed of traversal of the curve is unimportant) then one should consider the error measure
E(s, F) := itf { ],,F(T)
-s(n(i)),l’d~},
(19)
0
where (Yis an increasing homeomorphism We briefly consider the problem
of [0, 1 ] into itself.
inf{E (s, F) : Pi (s) E SF}
(20)
Under appropriate hypotheses of smoothness one can derive upper bounds on the error from the classical spline error bounds as follows. First notice that if Q is a diffeomorphism, we can make the change of variables 7 = Q (z) yielding
Thus if l/a’ < COwe can apply classical error estimates inf{(E(s,J’))
‘I2 : Pi(s) E S,“} < Kinf $ ,II(F a II I/
where K is a constant
depending
yielding 0 Q-‘)(~)JJ~ \tlk,
(22)
only on k.
5. Summary We hope that the ideas presented here provoke more discussion and experimentation. The general problem of recovering a geometric curve (or surface) when no explicit parameterization is given appears to be a pervasive problem. While we are unable to prove that ODR splines have favorable shape preserving properties, our numerical experiments suggest that they can play a very useful role in faithfully representing geometric shapes with fewer parameters than are generally required in competing methods.
S.P. Marin, P. W. Smith / Computer Aided Geometric Design 11 (1994) 247-267
261
Acknowledgement We would like to thank Weston Meyer for his insightful comments regarding uniqueness for ODR splines, and for suggesting an alternate, more direct counterexample to uniqueness.
References Boggs, P.T., Byrd, R.H. and Schnabel, R.B. (1987), A stable and efficient algorithm for nonlinear orthogonal distance regression, SIAM J. Sci. Stat. Comput. 8, 1052-1078. Boggs, P.T., Donaldson, J.T., Byrd, R.H. and Schnabel, R.B. (1989), ODRPACK: Software for weighted orthogonal distance regression, ACM Trans. Math. Software 15 (4), 348-364. de Boor, C. (1978), A Practical Guide to Splines, Springer, New York. Foley, T.A. and Nielson, GM. (1989), Knot selection for parametric spline interpolation, in Lyche, T. and Schumaker, L.L., eds., Mathematical Methods in Computer Aided Geometric Design, Academic Press, Boston, MA, 261-267. Herz, R.K. and Marin, S.P., (1980), Surface chemistry models of carbon monoxide oxidation on supported platinum catalysts, J. Catalysis 65, 281-296. Hoschek, J. (1988), Intrinsic parameterization for approximation, Computer Aided Geometric Design 5, 27-31. Lee, E.Y.T. ( 1989), Choosing nodes in parametric curve interpolation, Computer-Aided Design 21, 363-370. Marin, S.P. (1984), An approach to data parametrization in parametric cubic spline interpolation problems, J. Approx. Theory 41, 64-86. Meyer, W.W. (1978), Orthogonal projections. Solution to problem 988. Math. Magazine 51, 71-72. Meyer, W.W. (1990), Private communication. O’Dell, C.L. (1985), Approximating data with parametric B-splines for computer aided design, M.S. Thesis, Department of Computer Science, University of Utah, Salt Lake City. Pinkus, A. (1988), On the smoothest interpolation, SIAM J. Math. Anal. 19, 1431-1441. Plass, M. and Stone, M. (1983), Curve-fitting with piecewise parametric cubits, Computer Graphics 17, 229-239. Rademacher, C. and Scherer, K. (1990), Algorithms for computing best parametric cubic interpolation, in: Mason, J.C. and Cox, M.G., eds., Algorithms for Approximation II, Chapman and Hall, London, 193-207. Scherer, K. (1987), Uniqueness of best parametric interpolation by cubic spline curves, preprint. Scherer, K. and Smith, P.W. (1989), Existence of best parametric interpolation by curves, SIAM J. Math. Anal. 20, 160-168. Toepfer, H.J. (1981), Models for smooth curve fitting, Numererische Methoden der Approximationstheorie, 6, Internat. Ser. Numer. Math., Vol. 59, Birkhauser, Basel.