- Email: [email protected]
Surface fitting for interactive shape design
State-of-the-Art
Surface Fitting for Interactive Shape Design Joris S.M. Vergeest Faculty of Industrial Design Engineering, Delft University' of Technology, Jaffalaan 9, NL-2628 BX Delft, The Netherlands A method is presented to create smooth B-spline surfaces on basis of three-dimensional data points. The intended application is interactive shape design and shape representation with CAD workstations. This type of application allows specific assumptions to be made on the input data points, and it implies requirements in terms of speed and integration of the method with at least a graphics environment. There are no restrictions on the data locations, they need not fit into a grid and may be irregularly spaced. In an automatic process, the data are first represented by an array of B-spline curves. The parametrizations of these are harmonized as to enable approximate representation (within some tolerance) by a single B-spline surface with the smallest possible number of control points. The fitting system can also be used to represent a large amount of data points by an approximating surface, thus providing data compression.
Kevwords: CAD, Sculptured surfaces, B-splines, Shape design, Parametrization, Least squares fitting, Approximation, Computer graphics.
Joris S.M. Vergeest received his MS Degree in Physics and Computer Science in 1974 from the University of Nijmegen, The Netherlands. From the same institute he obtained his PhD in Experimental High Energy Physics in 1979. Dr. Vergeest is now scientist at Delft University of Technology, Faculty of Idustrial Design Engineering and is involved in research projects in the area of Geometric Modeling, CAD Data Exchange, Computer Graphics, Automated Robot Control and Automated NC. Elsevier Science Publishers B.V. Computers in Industry 13 (1989) 1-13 0166-3615/89/$3.50 © 1989 Elsevier Science Publishers B.V.
1. Introduction
Surface fitting techniques have been developed in a wide variety and for a great number of different applications. Interpolation and approximation schemes exist for three-dimensional scattered data as well as for data points with some predetermined structure, e.g. a rectangular grid arrangement. In addition, the fitting methods may vary with respect to the constraints which may be imposed, such as boundary conditions, smoothness criteria a n d / o r local deviation tolerances. Often, a surface fitting system is intended either for the representation of surfaces or for the design of shapes [2]. In the former case, the data to be fitted may be considered as stationary (e.g. measured geologic data) and the surface is mainly used to provide a data-compressed representation with a visually attractive form. In the case of shape design it may be expected that the data on which the fits are based will be interactively modified by a user. Depending on the requirements for the surface, interactive a n d / o r iterative shape manipulation techniques are used, e.g. on a CAD workstation. On B6zier/B-spline systems, operations like control point tweaking are getting possible with immediate display response. However, it is known that control point editing is not generally suited for fairing curves and surfaces. For B-spline surfaces, the tensor product concept hampers control point addition and deletion because a knot is associated with a row of the De Boor control polyhedron, rather than with a single vertex. Another difficulty is that judgement of surface quality involves ultimately computationally expensive presentations (e.g. Gaussian curvature). So, even when the geometric model is quickly updated, the graphical display takes a long time, slowing down the iterative process. On the other extreme, methods have been developed for non-interactive automatic surface construction, involving
2
State-of-the-Art
optimization constrained by fairness (expressed by some heuristic criterion), boundary derivatives and convexity preservation [1,9]. The fitting scheme presented here will be applicable to surface data representation as well as to design. The fitting facility has been embedded in an existing C A D / C A M system, SlPSURF [5], originating from Delft University and made commercially available by Delft Spline Systems, Delft, The Netherlands. SIPSURF has been used over the years for B-spline surface design but also to process and represent measured samples from various sources, with numbers of data points ranging up to 100,000. When the number of data points is very small (in the order of 100-200), surface fitting and display can be done in reasonably short time. Often, fast fitting of updated information is only feasible with local adaption methods. Our approach is based on a global method: change of one data point effects the surface as a whole. We will conclude that a global method is sufficiently fast for interactive design as long as the number of data points is small. In this article the surface fitting method is described. In the currently implemented version the surface is optimized with respect to deviation from the data points in a least-squares sense. In addition the surface will deviate from any data point less than some user-specified distance, while the number of control points is kept as small as possible. No additional criteria (e.g. overall smoothness) are taken into account by the fit. The idea is that such additional criteria should be checked interactively on the basis of high-quality images. Calculation times as measured for typical surface fits are reported, as well as graphical results obtained at the work station. Finally a list of improvements is given which are needed for either increasing the processing speed or for making the method suited for more general sets of data points.
2. SIPSURF and Tensor Product Surface Fitting
2.1. The SIPSURF System An a priori requirement for the fitting scheme is that the resulting surfaces are immediately usable in the SIPSURF modeler. This system is based on open and closed bivariate three-dimensional B-spline surfaces. The precise definition of the
C o m p u t e r s in I n d u s t r y
surfaces is as follows. In the three-dimensional Eucledian space, a B-spline surface S is defined by three coordinate functions of two common real variables (parameters) u and
S(u, ~,)= (x(u, ~), y(u, ~,), =(u, ,~,)), where x, y and z are bivariate functions: ,%,
lv~
/=1 .j-1
,%, x,
y(., .)= ?. Z t=l /=1 N,,
N,
i= 1 j=t
and (P~,,, Pv,, P = , ) = P i , are the three-dimensional control points of the surface. These control points are arranged in a matrix of size (N~ × N~,). N, and N~, are the number of control points in the u- and v-direction respectively. B,~o(u) is a normalized B-spline function of order k, (degree k , - 1 ) , defined on the knot vector for the u-parameter. This knot vector is a set of non-decreasing values u, ( i = 1, N , + k , ) , such that u,+k, > u, for i = 1, N,,. Then B,k,(u) can be calculated for u~, 4 u <~ ux, ~ ~ as follows: B]=
1, 0,
for u , ~ < u ~ < u , ~ , elsewhere;
and B,~o (u) -
u-u, bli+k,,
+
B,~"- ~(u
I - - Ut
Ui+k.
-- U
B~°i ~(u),
U i + k . - - lgt + 1
if k , > 1. The definition of B[,,(u) is completely analogous. Elaborate information on B-spline curves and surfaces can be found in [4] and references therein. SIPSURF employs specific knot vectors for open and closed surfaces. The four possible surface topologies are: (1) Surface open in u- and v-direction (planar topology).
Computers in Industry
J. Vergeest / Surface Fitting
(2) Surface open in u-direction and closed in o-direction (cylindrical topology). (3) Surface closed in u-direction and open in v-direction (cylindrical topology). (4) Surface closed in u- and v-direction (toroidal topology). If the surface is open in the u-direction, then the knot vector for the u-parameter is U1 . . . .
ua. + , = i ,
=Uk
~0
fori=l,
UN.+2 . . . . .
,
N.-k u+l,
UNu+ku = N. -- k . + 1.
This is a uniform knot sequence with knots of multiplicity k . at the parameter bounds. If the surface is closed in the u-direction, the knot distribution is strictly uniform: u,=i-k.,
i = 1 , N u + k ..
The same choice is made for the knot vector for the v-parameter. Further limits of the SIPSURF system are at present:
1 < N . , N,, < 100, 1 < k ~ , k~, < 22. In most cases SIPSURF is applied for shape design. The user has to determine somehow which control points P,j will produce the required surface S. Judgement of the shapes is done by means of computer-generated line drawings, color-shaded images or NC-milled prototypes [5,11]. As far as data fitting is concerned, two aspects of the specific choice of surfaces as described above are essential: (1) The B-spline surfaces are tensor product surfaces. On the one hand this gives a restriction to the types of data that can be adequately handled in a fit. On the other hand it offers an enormous reduction in terms of computation time and storage occupation compared to fitting of general surfaces. (2) The knot distributions are fixed and uniform. This means that optimization in some sense is excluded, as will be shown.
3
and reduced occupation of computer memory [8]. For the applications of shape design and surface representation the following can be assumed in a number of cases: (1) If the data is entered by a user for shape design, there will be points ordered in some way, such that contours or other characteristic (but otherwise arbitrary) curves are followed on or near the intended shape. (2) If the data originates from some measuring system, then the samples are taken in series (not necessarily along straight lines or planar curves). In either of these circumstances tensor product surface fitting is feasible. If the data is scattered without any derivable ordering, then the fitting method we describe will not be efficient~ as at least some gridding must be obtained first. In applications of industrial design however, irregularly scattered data typically do not occur. The principle of tensor product surface fitting is in its most simple form repeated curve fitting. Suppose that the data points are (or can be) arranged as M series of N points. Each series is approximated by a curve represented by K coefficients (control points). Then there are K series of M points available, each of which is approximated by a curve represented by L coefficients (using exactly the same fitting procedure as for the M series before). The resulting ( K × L) coefficients are exactly the coefficients (control points) of the intended approximating surface. Hence, rather than fitting one surface with ( K × L) variables to ( M × N ) data points, we now can split up the problem into ( M + K ) fits with K (or L) variables to N (or M ) data points. This offers an enormous efficiency gain in practice. Note that in our method it will be not required that the series of data points consist of the same number of points, nor that the data must be grouped into some grid. Also there is no restriction for the values taken by the coordinates.
3. Fitting B-Spline Surfaces to the Data--Description of the Method
2.2. Tensor Product Surface Fitting
If certain characteristics of the data are established, it is possible to benefit from the advantages of tensor product fitting: reduced calculation times
The principle of tensor product surface fitting described above will be worked out for the specific type of B-spline surface S ( u , v) depicted in Section 2.1.
4
State-of-the-Art
3.1. Assumptions Made for the Input Data There will be very few constraints on the starting information. The properties of the input data are: (1) The data consist of three-dimensional Cartesian points (i.e. values of x-, y- and z-coordinates). Two- or one-dimensional problems can be solved by setting one or two components to zero for all points. (2) The data points must be grouped in M series (2 ~< M ~< 99). The number 99 is arbitrary in as far as it represents a limit in the current software implementation. (3) Each series may consist of a different number ~ of data points D,j (i= l, M; j = l, N,) and 2 < ~ < 99 (same remark on the number 99 as above). (4) The order in which the data points occur in a series is relevant, as is the ordering of the series themselves. Normally the problem dictates the required data sequencing. On basis of these data a B-spline surface of order (k,, k,,) will be found that approximates each point Dsj within tolerance T = ( Tx, Ty, ~ ) in the sense that for each Dgj there exist a parameter value pair (u', v') such that
[S~(u', v')-D:~,j[<..Tx,
ls,,(u',
v') - :-,,,.l<- r,,,
.')If the N, are different for different i, the inequalities are not strictly true, as will be shown. In the current version of the method the order k, and k, are no variables in the fit, but have some fixed user-defined value smaller than 22 (see Section 2.1). T~, T~. and ~ may be assigned any non-negative value. Setting T = (0, 0, 0) would result in the interpolating B-spline surface that passes through all data points (exactly only if N, is independent of i).
3.2. Curve Fit to a Point Series The first stage of the method is a least-squares approximation of each point series by a B-spline curve r(u)= (x(u), y(u), z(u)) of order k , defined on a knot vector as given in Section 2.1. Consider one such point series consisting of N
('omputers in Industry'
data points D,, i = 1, N and N >~ k,. (Note that N >~k u is not a requirement for the fit; however in the case that 1 < N < k,, the point series problem is directed to a further stage of the scheme where underdetermined systems are solved; this will be described below.) Initially the x-, y- and :components of the D i are fitted independently by single-valued B-spline functions x(u), y(u) and z(u). The B-spline function x(u) is determined from the x-components of the D,: D~, . . . . . D~,. as follows. A B-spline function x ( u ) with the lowest possible number of control points (N,, = k,) is fitted to the data (the details of this fitting will be described below). If x(u) fails to approximate the values O~, within tolerance ½T~ then 5:, is increased by 1, and the fit is retried. If N, reaches the value N then the exact data interpolation will be obtained. Therefore a satisfactory fit with the smallest possible N, exists and has k u ~
x ( u ) = Y'~p~B,~"(u) i=l
to the values D~,, .... DxN proceeds as follows. The knot vector for the u-parameterization is set up as described in Section 2.1 for an open B-spline curve of order k, with N, control points. To each data point D,, a parameter value tj is assigned in the allowed u-range [uk, UN.+~]. This is necessary to formulate the problem as a fit of function x(u) to data points D,: such that N
2 w,(x(,,)- by,)2 /= 1
is minimal, where w: is a weight factor associated with data point D~. It should be remarked that in many approximation problems the t-values associated to the data points are known a priori, for which then some appropriate knot vector needs to be found [6], whereas in our scheme the data points and the knot vector are fixed, for which then some t-distribution needs to be found. The actual choice of the abscissae t/ is critical for the resulting x(u). A minimal requirement is that t/+~ > t / for I ~ < j < N and that there is at least one t in each non-zero interval (span)
Computers in Industry
J. Vergeest / Surface Fitting
[u i, u,+~] for k , < ~ i ~ N,. This still leaves many different possible t-distributions. If the total problem to be solved would be the fitting of curve x ( u ) or r(u), then strategies for optimal t-values could be followed. Similar optimizations are possible for data on regular grids. However, as we will see later, the t-distributions for the different point series must eventually be identical as must be the knot vectors of the resulting approximation curves. There are no spatial restrictions on the data points: one point series may be situated completely differently from any other one. This has led us to choose a rather regular spacing of t-values within the available knot intervals, taking into account dense or multiple knots. Suppose that a knot vector for the u-parameter has been set up for order k , and that the number of control points is N,. Then there are N, - k , + 1 knot intervals. At first, N,, locations s, are determined: U i + 1 ~'- b l i + 2 nk " " " 4- U i + k , ,
si=
l
k.-1
i=1,
N.,.
The s, are scaled and shifted to produce s,' such that s 1 and s,v are at the interval extremes u~. and UN.+~" t
s'=sJR-
D,
will be distributed such that there exists at least one in each non-zero interval, of which there are N u - k , + 1 (or less if multiple interior knots would be allowed). The weights % in (1) are all taken identically 1; they play no role in the scheme described here. In Section 5 a possible usage of data point weights for tensor product surface fitting will be indicated. With the knot vector and t-distribution as described above, it can be verified that there exists a unique set P,,, . . . . P,, that solves (1). With all w.I= 1 this means that"{he vector
B("(t,
...
B~::(t,)
P,,
D,,
...
B~,':( t u )
P~,
D,,
s o(t ) B~" (t v )
has minimal length (minimal sum of squares of the components). The matrix of size (N,, × N ) above contains many zero's if k , << N,, << N, due to the local support of the B-spline functions:
t
i= l, N,,
where R
5
B f " ( o ) = O ,
i f t , ~ [ u , ,
ui+k. I.
As a result of the ordering of the matrix rows according to increasing t, the non-zero elements are grouped in a sense depicted below.
=
X X X X
D = s l / R - uk.,. The number of t-values ( N ) is larger than or equal to N.; therefore if
A =(N,-1)/(N-1), then 0 < A ~< 1. Finally, t~ is taken as follows:
. . . . . . . . . .
- X X X X
. . . . . . . . . .
• . X X X X
.........
• . X X X X
.........
• . . X X X X
1) + 1 - r n ) ( s ' + , - s2) ,
........
....
X X X X
.......
....
X X X X
.......
.....
tj = s[. + ( A ( j -
. . . . . . . . . . .
• X X X X
......
X X X X X X X X
...... .....
where m = INT(A(j-
1)) + 1
and I N T ( a ) is the largest integer less than or equal to a. This guarantees that t 1 and t~, are located at the interval bounds, which is supposedly an appropriate condition for the fit. Further, if N, approaches N, the t-distribution will approach the s-distribution as obtained by so-called knot averaging, which is recommended for interpolation (the case N, = N ) [6]. In any case the N t-values
where x stands for some positive value and • for zero. The (N, X N , ) coefficient matrix for the normal equations to solve (1) turns out to be positive semi-definite, symmetric and to have a band width less than k , [6]. This makes the system suited for solution by Choleski factorization using De Boor's routines L2APPR, BCHFAC and BCHSLV. The same procedure is followed to obtain y ( u ) and z(u) to approximate the D,. and D. within tolerance ~ ~ and ½T: respective(y. In general this would lead to different numbers of con]
•
~;
6
State-of-the-Art
Computers in Indust~
ill
i
o
i
l
1
I ill 2
3
1tll I
J~J~____J_AL
0
2
!
3
4
S
6
7
: ...... :7
/ 3
fD .
__
!
[]
t
,
'il
/ /
[]
data points
•
curve points
L-]i'
j,
\
/
T× = lo i
L. . . . .
I ......
0
LENGHT
50
~,
=
10
d
100
UNITS
jJ :),
jl /~ //
k
¢
' i
~[
:i ~ /
/
L
i'
i
e
k
i
d
L
Fig. 1. Refitting a B-spline curve on a new uniform knot vector with a larger number of knots, while avoiding to get closer to the data points than necessary. Example projected onto the XY-plane. (a) 15 data points, B-spline curve, its control polygon and its knot vector to approximate the data within tolerance 7~, = 10 and T~. = 10 length units, with minimal number of control points (N~ = 6). Suppose that fits for the other curves imply that the number of control points must be increased to 10; then 10 abscissae are determined in the parameter range as described in the text. These and the corresponding curve points have been calculated and are shown in the figure. (b) B-spline curve with its control polygon that interpolates the 10 points evaluated on the curve of Fig. l(a). In the new knot vector the 10 abscissae used for interpolation are indicated. (c) Comparison between the curve of Fig. l(a) (indicated by arrows) and the refit of Fig. l(b). (d) Comparison between the curve of Fig. l(a) (indicated by arrows) and the B-spline curve that would be found by refitting the original data with 10 control points.
trol points for the x-, y- and z - c o m p o n e n t s , which w o u l d not produce a curve r(u) with three-dimensional control points. Therefore, in the i m p l e m e n tation the tolerance checks are carried out for all
three c o m p o n e n t s , for each N . . If o n e or m o r e of the c o m p o n e n t s deviate too m u c h from the data, N. is incremented for all three fits with x ( u ) , y(u) and z(u).
Cornputers in lnd~t~
J. Vergeest / Surface Fitting
At this point we have obtained a B-spline curve of order k,: U, = E i=1
that approximates ones row of data points in the least-squares sense, and in addition does not deviate from the data by more than half of the specified tolerance T.
3.3. Re-Representation of a B-Spline Curve on New Uniform Knots For each of the M rows of data points a B-spline curve ~)(u) can be determined in the way described above: N,,j
rj(u) = Z P, j B [ " ( u ) ,
j = 1 . . . . . M.
Note that the points P,j will in general not form a rectangular array. To stay within tolerance ±2 7", a different number N,, of control points may be needed for the different rows, depending e.g. on the arbitrary number of points contained in the rows. On the other hand, in the following stage of tensor product surface fit, we will need curves rj(u) on identical knot vectors and hence with identical N,; The obvious choice is to take N,m., = maximum of N~,, .... N,M, which implies that each N,, must be increased by (N, .... - N , , ) . There are several ways to do this. One way is to replace N,, by N,,,,. in the curve fit for r(u). However, if N,m.. > Nt for some j ( N / i s the number of data points in row j), the systems gets underdetermined. Even if N, .... ~< Nj, the newly fitted curve r(u) would have the disadvantage of being less smooth than necessary for the required tolerance (see Fig. l(c)). Therefore, we have chosen to keep the geometry of the curves r~(u) principally unchanged while increasing N , . This is accomplished in the following way. N,,,.... locations e,j on curve ~(u) are evaluated at parameter values t,, i = 1. . . . . N,m., e,, = ~/(t,),
where t, is determined in the parameter range of r(u) as described above, except that N (the number of abscissae) is replaced by N.m.. Then the locations e,j are regarded as N~..... data points to
7
which a B-spline curve of order k , with N u...... control points has to be fitted (which is in this case equivalent to interpolation). This too is performed in the way outlined above: the appropriate uniform knot vector is set up and Num~. abscissae are placed in the parameter range. The fitted curve will exactly interpolate the N, ..... data points and it can be verified that under the actual conditions of knots and t-values, the deviation of the curve from the point set r~(u) is small compared to the tolerance T, except for some extreme cases with large C k,,- ~ discontinuities. A detailed analysis would be required on to how ~(u) can be represented optimally under the given constraints. Note that increasing the number of control points from N~j to N. cannot be simply performed by knot insertion. This would in general produce a non-uniform knot vector, which is not allowed by the embedding SIPSURF system. Our re-representation is on basis of a new, but still uniform knot vector, for which in general an exact solution does not exist. Any underdetermined system (curve with N/< k . ) possibly encountered earlier in the scheme, will be solved here. If ~ < k . then the exactly interpolating B-spline curve ~(u) of order Nj is determined, on which N...~ locations are evaluated in the way just described, to which ~(u) of order k. with N.m.. control points is fitted.
3.4. Obtaining a Surface from the Fitted Curves We have now obtained M B-spline curves r~(u) of order k. with N.m., control points E,I each ( i--- 1 . . . . . Nu.,~ ; j = 1 . . . . . M ). Next we regard these control points as N,....... rows each consisting of M points. In principle the procedures described in the previous two sections are repeated, but now in the inverse direction. For every row i we determine by least-squares approximation the Bspline curve r,(v) of order k, that approaches the M points E,/ ( j = 1. M ) within tolerance ½T with the minimal number of control points N,,. By using exactly the same knot and abscissae conventions as for the curve fits described above and by reparametrizations, we will obtain N, ..... B-spline curves ~(v) with identical numbers (N,.m.) of control points P,j (i = 1 . . . . . N ...... ; j = 1 . . . . . N,,m.,). The number N , , will depend on the data, but it will be less than or equal to M and is the largest N,, encountered. Due to properties of
8
State-of-the-Art
Uomputers in Industry
Fig. 2. Surface fits with tolerance T = (0, 0, 0). (a) Data points and B-spline surface fit with degree = 1 (k,, = k,. = 2); (b) degree = 2; (c) degree = 3; (d) degree ~ 9; (e) overshoot in the fitted surface with degree = 9.
Computers in Industry
J. Vergeest / Surface Fitting
Table 1 Fits of a bicubic B-spline surface ( k , = k,, = 4) to three different p o i n t sets at different tolerances N u m b e r of d a t a rows
N u m b e r of data points
Tolerance
N u m b e r of coefficients
6 6 5 5 5 35 35 35 35
21 21 160 160 160 1775 1775 1775 1775
1.0 10.0 1.0 10.0 20.0 1.0 10.0 20.0 100.0
5X 6 = 5× 6 = 36x 5 = 28 X 5 = 20x 4 = 61 × 35 = 43 × 34 = 32 × 34 = 4×20 =
30 30 180 140 80 2135 1462 1088 80
Elapsed time a (seconds) SFIT01
(general)
1 1 10 7 4 293 150 87 12
(15) (15) (8870) (1200) (350) (-) (-) (23600) (955)
Elapsed times are as m e a s u r e d for runs on a Digital VAX-Station 3200.
tensor product surfaces, the control points P,/ are the (N,m,~ × N,,m,~) control points for the B-spline surface S(u, v) of order ( k , , k,,) defined on the former knot vectors for the u- and v-parameter that approximates the points E u in the leastsquares sense. As the curves ~)(u) deviate at most approximately ½T from the original data points D,j, while the curves r,(v) are closer than ~ T to the control points of the curves ~)(u), the surface S(u, v) approximates the original data within T, except for perhaps some extreme cases. In the current implementation there is no explicit distance check between D,j and S.
4. Implementation The fitting system has been implemented in standard Fortran as a computer program named SFIT0I. Interfacing with the interactive SIPSURF system proceeds, momentarily, through data base level operations. The user produces, somehow, a list of coordinate values of the data points in free format. The data may originate from a 2-D or 3-D digitizer, or be entered directly into a text file. The desired surface degree in u- and v-direction and the x-, y- and z-tolerance must be specified. Following the method described above, SFIT01 adds the new B-spline surface to the SIPSURF database. The user gets a report of which data points deviate maximally from the surface, and in which direction. It is the intention to arrive at a full integration of SFIT01 with the workstation graphics in the future.
4.1. Processing times" In the current implementation SFIT0I spends most of its time on finding the smallest knot sets, rather than on calculating the final control points. The program can be significantly optimized in the knot searching part. Still, as long as the number of data points is small, the total time to find the surface that passes the tolerance tests is short compared to the time needed for raster image generation on the particular work station. Table 1 lists the elapsed time in seconds to produce a SIPSURF B-spline surface from the data points. The degree in both directions is set to 3 in all cases. The time needed by SFIT01 (which uses De Boor's L 2 A P P R ) has been compared to the time spent by a similar program in which no account is taken of the specific local support properties of B-spline curves, but where instead a general least-squares solver is applied (routine E04FDF from N A G Library [10]). The programs have been run on a Digital's VAX-Station 3200 including a graphical display of 1024 × 768 pixels with 8 bit planes. Generation of a high-lighted color picture in a 256 × 256 pixel window typically takes 60 seconds. In the next section the possibilities of integration of SFIT01 with SIPSURF to enable interactive fitting will be indicated. As can be seen in Table 1, the number of coefficients found for the B-spline surface can be larger than the total number of data points when the tolerance is small. This is due to the fact that a row consisting of few data points is eventually fitted with a curve having as many control points
10
State-of-the-Art
('omputers in Industry
Fig. 3. Fits with B-spline degree = 2 to the data points of Fig, 2 for different values of the tolerance T: (a) 100 cm: (b) 50 cm: (c) 30 cm; (d) 20 cm; (e) 10 cm; (f) 1 cm.
Computers in Industry
a
J. Vergeest / Surface Fitting
11
b
Fig. 4. Highlighted image (left) and NC milled part (right; scale 1 : 20) from the model of Fig. 2(c). The symmetric other half of the car has been added to the model by using SlPSURVoperations.
as the maximum needed for any other data row (Section 3.3).
4.2. Example As an example, a set of very few data points (a total of 75 divided over 20 rows) has been taken from a simple four-view drawing of a car model (advertisement leaflet from the manufacturer). In each row only 3 to 5 data points were determined on the left half of the car body. Figure 2 shows the data points and fits to the data at zero-tolerance with B-spline degree 1, 2, 3 and 9. For degree 1, 2 and 3, B-spline surfaces with 20 × 5 control points were obtained, for degree 9 one with 20 × 10 control points. The gross shape of the car body is reproduced with one single B-spline surface. The fit with degree 2 seems optimal, as it features the fewest overshoots (curling self-intersecting surface). To remove such overshoots, the user may add a few data points to guide the surface fit. Alternatively, the tolerance can be increased. In Fig. 3 the results of fits with degree 2 for different tolerance values are shown. B-spline surfaces of 20 × 5 control points are found for tolerances up to 20 cm (length of car is approximately 400 cm). For larger tolerances fewer control points are sufficient. Figure 4 shows an NC milled model out of foam material (after adding symmetrically the other half of the car body using SlPSUaF oper-
ations " G l u e " a n d / o r "Connec" for automatic G 2 blending). The question arises of what the user should do to further modify and detail the shape. A standard option of SlPSUaF is control point editing (designation with a mouse and typing delta coordinate values). However, the control points found by fits are not easily interpretable generally (Fig. 5). It will be investigated in the future what the possibilities are of moving and adding data points (rather than control points) in an interactive way.
5. Conclusions and Further Research
5.1. Conclusions A method has been presented to generate a single B-spline surface on the basis of (possibly very few) data points. Smoothness and closeness to the data is controlled by simple tolerance values. The program SliT01 automatically finds the B-spline surface with the smallest number of control points to meet the tolerance requirement. We have illustrated that it would not have been possible to find those control points in a classical trial-and-error procedure, directly by the user. This implies that further control point editing after the fit is not practical in general. On the other hand, since the fitting procedure is sufficiently fast for
12
State-of-the-Art
Computers in lndustrv
-: j )
Fig. 5. The control points found by SF~T01are often not suited for direct further modeling by the user (shown are the control points of the surface of Fig. 2(c)).
less than about 300 data points, interactive modifications, addition and deletion of single Or groups of data points seem feasible in implementations on present workstations. 5.2. Further Points o f Interest
There are alternatives to some of the choices made for the fitting procedure. Moreover, several improvements and extensions to the currently implemented system are possible. Some of the issues are: (1) Most of the computation time in SFIT01 is spent on finding the simplest possible B-spline surface, not on fitting its control points. Much time could be safed if, instead of simply incrementing N, by one (Section 3.2) a more direct search strategy (e.g. a binary search) would be followed to obtain the optimal N, for the curve fits. (2) The maximal deviation of the surface from the data points is the criterion for SFIT01 to accept the fit or to increase the number of B-spline control points. A direct way to check against this criterion is possibly with surfaces that minimize the L ~ norm (largest deviation). Such fits can be implemented with N A G Library routine E02GCF, but these are as slow as the general solver (E04FDF, see Section 4.1). We could not find an L ~ minimizer for B-spline curves that could approach the speed of the De Boor's L 2 (least-squares) algorithm, as implemented in SFIT01. (3) SFIT01 can be extended to closed and doubly closed B-spline surfaces (the three other topologies of SIPSURF, see Section 2.1). This would be
very appropriate if the data points are known to be located on closed contours. For closed (periodic) curves, the matrix of size iN,, × N ) in Section 3.2 will contain a few more non-zero elements in the top-right a n d / o r in the left-bottom corners, causing the loss of the band structure of the normal equations matrix. The system can still be solved with the general N A G Library routine E04FDF. To take advantage of the (still highly) structured matrix, algorithms due to Dierckx can be applied [7]. (4) The parametrization of any fitted curve is independent from the geometry of nearby other curves (except the size of the knot vector). Suppose that the M curves found in Section 3.4 happen to be well alligned or even parallel in a plane. If one of those curves was obtained from a row of data points with a large concentration of points at a particular place, then there will be also a concentration of control points of the fitted B-spline curve near that place. If nearby data rows do not exhibit such a concentration, the surface as a whole, found in Section 3.4 may be distorted. For some applications this may just be the desired effect. In many cases, the surface should depend on the geometry of the intermediate curves, not on their parametrizations. Parametrization independence can be obtained by evaluating points on a curve at intervals of equal curve length and then refitting the curve to those points. This could be an alternative for the re-representation described in Section 3.3. For fits to data on closed contours, the reparametrization should even allow shifts of
Fig. 6. Unwanted twists and folds in the surface if the parametrization does not account for irregularly spaced data points.
J. Vergeest / Surface Fitting
Computers in Industry
the begin/end locations of the curves to avoid unwanted twists in the resulting surface. Figure 6 gives an example of a surface obtained by SFIT0I from uncorrected photogrammetric contour data. (5) The curve fits in SFIT01 are unweighted. If there is knowledge about relative weights of the data points, it is straightforward to set of w s in (1) accordingly. It is not clear what the effect of the weights should be in the further stages of the procedure. One possibility is to propagate the data point weights to the relevant curve control points for the transverse fits. (6) The modeling capacity of SFIT0I when integrated with the SIPSURr surface modeler needs further investigation. When raster image display is sufficiently fast, will then the user be able to obtain shapes with all necessary features? Unwanted fluctuations should be damped by adding or removing data points, or by changing their weights. (7) SFIT0I can be significantly speeded up for operations that effect only one data point or one row of data points (these are the operations most likely to occur interactively). The curve fits on the other rows then need not be redone in many cases. The transverse fits must all be repeated, but a large part of the matrix equations remains unchanged. (8) User interface and display graphics should be optimized for interactive shape design. Matters of interest are movement of data points in 3-D, fast coarse image generation with automatic refinement until the next user interrupt, storage of a number of recent images including the data from which the user operations causing those images can be deduced.
Acknowledgements The author wishes to thank Wim Sybranda for making available preliminary design data, Rinus
13
Grabijn for preparing illustrations and Bram de Smit for the production of the NC-model.
References [1] E. Andersson, R. Andersson, M. Boman, T. Elmroth, B. Dahlberg and B. Johansson, "Automatic construction of surfaces with prescribed shape", Comput. Aided Des., Vol. 20, No. 6, 1988. [2] R.E. Barnhilt, "Representation and approximation of surfaces", in: J.R. Rice, ed., Mathematical Software 111, Academic Press, New York, 1977, pp. 69-120.' [3] R.E. Barnhill, "'A survey of the representation and design of surfaces", 1EEE Comput. Graphics AppL, Vol. 3, No. 10, October 1983, pp. 9-16. [4] W. Boehm, G. Farin and J. Kahmann, "A survey of curve and surface methods in CAGD", Comput. Aided Geom. Des., Vol. 1, No. 1, 1984, pp. 1-60. [5] J.J. Brock, "CADAMP, Computer Aided Design and Model Production", Proc. DECUS '84, Digital Equipment Corporation, Marlboro, MA, 1984. [6] C. de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978. [7] P. Dierckx, FITPACK User Guide, Part 1: Curve Fitting Routines, Report TW 89, Dept. of Computer Science, Katholieke Universiteit Leuven, Belgium, 1987. [8] J.G. Hayes, "Fitting surfaces to data", in: R.R. Martin, ed., The Mathematics of Surfaces 11, Clarendon Press, Oxford, 1987, pp. 17-38. [9] N.J. Lott and D.1. Pullin, "Method for fairing B-spline surfaces", Comput. Aided Des., Vol. 20, Nr. 10, 1988, pp. 597-604. [10] NAG, The NAG Fortran Library Manual--Mark 12. Numerical Algorithms Group, Oxford, 1987. [11] J.S.M. Vergeest, "Connecting arbitrary surfaces under geometric constraints", Computers in lndustlT, Vol. 8, No. 1, 1987, pp. 3-12.
Surface Fitting for Interactive Shape Design Joris S.M. Vergeest Faculty of Industrial Design Engineering, Delft University' of Technology, Jaffalaan 9, NL-2628 BX Delft, The Netherlands A method is presented to create smooth B-spline surfaces on basis of three-dimensional data points. The intended application is interactive shape design and shape representation with CAD workstations. This type of application allows specific assumptions to be made on the input data points, and it implies requirements in terms of speed and integration of the method with at least a graphics environment. There are no restrictions on the data locations, they need not fit into a grid and may be irregularly spaced. In an automatic process, the data are first represented by an array of B-spline curves. The parametrizations of these are harmonized as to enable approximate representation (within some tolerance) by a single B-spline surface with the smallest possible number of control points. The fitting system can also be used to represent a large amount of data points by an approximating surface, thus providing data compression.
Kevwords: CAD, Sculptured surfaces, B-splines, Shape design, Parametrization, Least squares fitting, Approximation, Computer graphics.
Joris S.M. Vergeest received his MS Degree in Physics and Computer Science in 1974 from the University of Nijmegen, The Netherlands. From the same institute he obtained his PhD in Experimental High Energy Physics in 1979. Dr. Vergeest is now scientist at Delft University of Technology, Faculty of Idustrial Design Engineering and is involved in research projects in the area of Geometric Modeling, CAD Data Exchange, Computer Graphics, Automated Robot Control and Automated NC. Elsevier Science Publishers B.V. Computers in Industry 13 (1989) 1-13 0166-3615/89/$3.50 © 1989 Elsevier Science Publishers B.V.
1. Introduction
Surface fitting techniques have been developed in a wide variety and for a great number of different applications. Interpolation and approximation schemes exist for three-dimensional scattered data as well as for data points with some predetermined structure, e.g. a rectangular grid arrangement. In addition, the fitting methods may vary with respect to the constraints which may be imposed, such as boundary conditions, smoothness criteria a n d / o r local deviation tolerances. Often, a surface fitting system is intended either for the representation of surfaces or for the design of shapes [2]. In the former case, the data to be fitted may be considered as stationary (e.g. measured geologic data) and the surface is mainly used to provide a data-compressed representation with a visually attractive form. In the case of shape design it may be expected that the data on which the fits are based will be interactively modified by a user. Depending on the requirements for the surface, interactive a n d / o r iterative shape manipulation techniques are used, e.g. on a CAD workstation. On B6zier/B-spline systems, operations like control point tweaking are getting possible with immediate display response. However, it is known that control point editing is not generally suited for fairing curves and surfaces. For B-spline surfaces, the tensor product concept hampers control point addition and deletion because a knot is associated with a row of the De Boor control polyhedron, rather than with a single vertex. Another difficulty is that judgement of surface quality involves ultimately computationally expensive presentations (e.g. Gaussian curvature). So, even when the geometric model is quickly updated, the graphical display takes a long time, slowing down the iterative process. On the other extreme, methods have been developed for non-interactive automatic surface construction, involving
2
State-of-the-Art
optimization constrained by fairness (expressed by some heuristic criterion), boundary derivatives and convexity preservation [1,9]. The fitting scheme presented here will be applicable to surface data representation as well as to design. The fitting facility has been embedded in an existing C A D / C A M system, SlPSURF [5], originating from Delft University and made commercially available by Delft Spline Systems, Delft, The Netherlands. SIPSURF has been used over the years for B-spline surface design but also to process and represent measured samples from various sources, with numbers of data points ranging up to 100,000. When the number of data points is very small (in the order of 100-200), surface fitting and display can be done in reasonably short time. Often, fast fitting of updated information is only feasible with local adaption methods. Our approach is based on a global method: change of one data point effects the surface as a whole. We will conclude that a global method is sufficiently fast for interactive design as long as the number of data points is small. In this article the surface fitting method is described. In the currently implemented version the surface is optimized with respect to deviation from the data points in a least-squares sense. In addition the surface will deviate from any data point less than some user-specified distance, while the number of control points is kept as small as possible. No additional criteria (e.g. overall smoothness) are taken into account by the fit. The idea is that such additional criteria should be checked interactively on the basis of high-quality images. Calculation times as measured for typical surface fits are reported, as well as graphical results obtained at the work station. Finally a list of improvements is given which are needed for either increasing the processing speed or for making the method suited for more general sets of data points.
2. SIPSURF and Tensor Product Surface Fitting
2.1. The SIPSURF System An a priori requirement for the fitting scheme is that the resulting surfaces are immediately usable in the SIPSURF modeler. This system is based on open and closed bivariate three-dimensional B-spline surfaces. The precise definition of the
C o m p u t e r s in I n d u s t r y
surfaces is as follows. In the three-dimensional Eucledian space, a B-spline surface S is defined by three coordinate functions of two common real variables (parameters) u and
S(u, ~,)= (x(u, ~), y(u, ~,), =(u, ,~,)), where x, y and z are bivariate functions: ,%,
lv~
/=1 .j-1
,%, x,
y(., .)= ?. Z t=l /=1 N,,
N,
i= 1 j=t
and (P~,,, Pv,, P = , ) = P i , are the three-dimensional control points of the surface. These control points are arranged in a matrix of size (N~ × N~,). N, and N~, are the number of control points in the u- and v-direction respectively. B,~o(u) is a normalized B-spline function of order k, (degree k , - 1 ) , defined on the knot vector for the u-parameter. This knot vector is a set of non-decreasing values u, ( i = 1, N , + k , ) , such that u,+k, > u, for i = 1, N,,. Then B,k,(u) can be calculated for u~, 4 u <~ ux, ~ ~ as follows: B]=
1, 0,
for u , ~ < u ~ < u , ~ , elsewhere;
and B,~o (u) -
u-u, bli+k,,
+
B,~"- ~(u
I - - Ut
Ui+k.
-- U
B~°i ~(u),
U i + k . - - lgt + 1
if k , > 1. The definition of B[,,(u) is completely analogous. Elaborate information on B-spline curves and surfaces can be found in [4] and references therein. SIPSURF employs specific knot vectors for open and closed surfaces. The four possible surface topologies are: (1) Surface open in u- and v-direction (planar topology).
Computers in Industry
J. Vergeest / Surface Fitting
(2) Surface open in u-direction and closed in o-direction (cylindrical topology). (3) Surface closed in u-direction and open in v-direction (cylindrical topology). (4) Surface closed in u- and v-direction (toroidal topology). If the surface is open in the u-direction, then the knot vector for the u-parameter is U1 . . . .
ua. + , = i ,
=Uk
~0
fori=l,
UN.+2 . . . . .
,
N.-k u+l,
UNu+ku = N. -- k . + 1.
This is a uniform knot sequence with knots of multiplicity k . at the parameter bounds. If the surface is closed in the u-direction, the knot distribution is strictly uniform: u,=i-k.,
i = 1 , N u + k ..
The same choice is made for the knot vector for the v-parameter. Further limits of the SIPSURF system are at present:
1 < N . , N,, < 100, 1 < k ~ , k~, < 22. In most cases SIPSURF is applied for shape design. The user has to determine somehow which control points P,j will produce the required surface S. Judgement of the shapes is done by means of computer-generated line drawings, color-shaded images or NC-milled prototypes [5,11]. As far as data fitting is concerned, two aspects of the specific choice of surfaces as described above are essential: (1) The B-spline surfaces are tensor product surfaces. On the one hand this gives a restriction to the types of data that can be adequately handled in a fit. On the other hand it offers an enormous reduction in terms of computation time and storage occupation compared to fitting of general surfaces. (2) The knot distributions are fixed and uniform. This means that optimization in some sense is excluded, as will be shown.
3
and reduced occupation of computer memory [8]. For the applications of shape design and surface representation the following can be assumed in a number of cases: (1) If the data is entered by a user for shape design, there will be points ordered in some way, such that contours or other characteristic (but otherwise arbitrary) curves are followed on or near the intended shape. (2) If the data originates from some measuring system, then the samples are taken in series (not necessarily along straight lines or planar curves). In either of these circumstances tensor product surface fitting is feasible. If the data is scattered without any derivable ordering, then the fitting method we describe will not be efficient~ as at least some gridding must be obtained first. In applications of industrial design however, irregularly scattered data typically do not occur. The principle of tensor product surface fitting is in its most simple form repeated curve fitting. Suppose that the data points are (or can be) arranged as M series of N points. Each series is approximated by a curve represented by K coefficients (control points). Then there are K series of M points available, each of which is approximated by a curve represented by L coefficients (using exactly the same fitting procedure as for the M series before). The resulting ( K × L) coefficients are exactly the coefficients (control points) of the intended approximating surface. Hence, rather than fitting one surface with ( K × L) variables to ( M × N ) data points, we now can split up the problem into ( M + K ) fits with K (or L) variables to N (or M ) data points. This offers an enormous efficiency gain in practice. Note that in our method it will be not required that the series of data points consist of the same number of points, nor that the data must be grouped into some grid. Also there is no restriction for the values taken by the coordinates.
3. Fitting B-Spline Surfaces to the Data--Description of the Method
2.2. Tensor Product Surface Fitting
If certain characteristics of the data are established, it is possible to benefit from the advantages of tensor product fitting: reduced calculation times
The principle of tensor product surface fitting described above will be worked out for the specific type of B-spline surface S ( u , v) depicted in Section 2.1.
4
State-of-the-Art
3.1. Assumptions Made for the Input Data There will be very few constraints on the starting information. The properties of the input data are: (1) The data consist of three-dimensional Cartesian points (i.e. values of x-, y- and z-coordinates). Two- or one-dimensional problems can be solved by setting one or two components to zero for all points. (2) The data points must be grouped in M series (2 ~< M ~< 99). The number 99 is arbitrary in as far as it represents a limit in the current software implementation. (3) Each series may consist of a different number ~ of data points D,j (i= l, M; j = l, N,) and 2 < ~ < 99 (same remark on the number 99 as above). (4) The order in which the data points occur in a series is relevant, as is the ordering of the series themselves. Normally the problem dictates the required data sequencing. On basis of these data a B-spline surface of order (k,, k,,) will be found that approximates each point Dsj within tolerance T = ( Tx, Ty, ~ ) in the sense that for each Dgj there exist a parameter value pair (u', v') such that
[S~(u', v')-D:~,j[<..Tx,
ls,,(u',
v') - :-,,,.l<- r,,,
.')If the N, are different for different i, the inequalities are not strictly true, as will be shown. In the current version of the method the order k, and k, are no variables in the fit, but have some fixed user-defined value smaller than 22 (see Section 2.1). T~, T~. and ~ may be assigned any non-negative value. Setting T = (0, 0, 0) would result in the interpolating B-spline surface that passes through all data points (exactly only if N, is independent of i).
3.2. Curve Fit to a Point Series The first stage of the method is a least-squares approximation of each point series by a B-spline curve r(u)= (x(u), y(u), z(u)) of order k , defined on a knot vector as given in Section 2.1. Consider one such point series consisting of N
('omputers in Industry'
data points D,, i = 1, N and N >~ k,. (Note that N >~k u is not a requirement for the fit; however in the case that 1 < N < k,, the point series problem is directed to a further stage of the scheme where underdetermined systems are solved; this will be described below.) Initially the x-, y- and :components of the D i are fitted independently by single-valued B-spline functions x(u), y(u) and z(u). The B-spline function x(u) is determined from the x-components of the D,: D~, . . . . . D~,. as follows. A B-spline function x ( u ) with the lowest possible number of control points (N,, = k,) is fitted to the data (the details of this fitting will be described below). If x(u) fails to approximate the values O~, within tolerance ½T~ then 5:, is increased by 1, and the fit is retried. If N, reaches the value N then the exact data interpolation will be obtained. Therefore a satisfactory fit with the smallest possible N, exists and has k u ~
x ( u ) = Y'~p~B,~"(u) i=l
to the values D~,, .... DxN proceeds as follows. The knot vector for the u-parameterization is set up as described in Section 2.1 for an open B-spline curve of order k, with N, control points. To each data point D,, a parameter value tj is assigned in the allowed u-range [uk, UN.+~]. This is necessary to formulate the problem as a fit of function x(u) to data points D,: such that N
2 w,(x(,,)- by,)2 /= 1
is minimal, where w: is a weight factor associated with data point D~. It should be remarked that in many approximation problems the t-values associated to the data points are known a priori, for which then some appropriate knot vector needs to be found [6], whereas in our scheme the data points and the knot vector are fixed, for which then some t-distribution needs to be found. The actual choice of the abscissae t/ is critical for the resulting x(u). A minimal requirement is that t/+~ > t / for I ~ < j < N and that there is at least one t in each non-zero interval (span)
Computers in Industry
J. Vergeest / Surface Fitting
[u i, u,+~] for k , < ~ i ~ N,. This still leaves many different possible t-distributions. If the total problem to be solved would be the fitting of curve x ( u ) or r(u), then strategies for optimal t-values could be followed. Similar optimizations are possible for data on regular grids. However, as we will see later, the t-distributions for the different point series must eventually be identical as must be the knot vectors of the resulting approximation curves. There are no spatial restrictions on the data points: one point series may be situated completely differently from any other one. This has led us to choose a rather regular spacing of t-values within the available knot intervals, taking into account dense or multiple knots. Suppose that a knot vector for the u-parameter has been set up for order k , and that the number of control points is N,. Then there are N, - k , + 1 knot intervals. At first, N,, locations s, are determined: U i + 1 ~'- b l i + 2 nk " " " 4- U i + k , ,
si=
l
k.-1
i=1,
N.,.
The s, are scaled and shifted to produce s,' such that s 1 and s,v are at the interval extremes u~. and UN.+~" t
s'=sJR-
D,
will be distributed such that there exists at least one in each non-zero interval, of which there are N u - k , + 1 (or less if multiple interior knots would be allowed). The weights % in (1) are all taken identically 1; they play no role in the scheme described here. In Section 5 a possible usage of data point weights for tensor product surface fitting will be indicated. With the knot vector and t-distribution as described above, it can be verified that there exists a unique set P,,, . . . . P,, that solves (1). With all w.I= 1 this means that"{he vector
B("(t,
...
B~::(t,)
P,,
D,,
...
B~,':( t u )
P~,
D,,
s o(t ) B~" (t v )
has minimal length (minimal sum of squares of the components). The matrix of size (N,, × N ) above contains many zero's if k , << N,, << N, due to the local support of the B-spline functions:
t
i= l, N,,
where R
5
B f " ( o ) = O ,
i f t , ~ [ u , ,
ui+k. I.
As a result of the ordering of the matrix rows according to increasing t, the non-zero elements are grouped in a sense depicted below.
=
X X X X
D = s l / R - uk.,. The number of t-values ( N ) is larger than or equal to N.; therefore if
A =(N,-1)/(N-1), then 0 < A ~< 1. Finally, t~ is taken as follows:
. . . . . . . . . .
- X X X X
. . . . . . . . . .
• . X X X X
.........
• . X X X X
.........
• . . X X X X
1) + 1 - r n ) ( s ' + , - s2) ,
........
....
X X X X
.......
....
X X X X
.......
.....
tj = s[. + ( A ( j -
. . . . . . . . . . .
• X X X X
......
X X X X X X X X
...... .....
where m = INT(A(j-
1)) + 1
and I N T ( a ) is the largest integer less than or equal to a. This guarantees that t 1 and t~, are located at the interval bounds, which is supposedly an appropriate condition for the fit. Further, if N, approaches N, the t-distribution will approach the s-distribution as obtained by so-called knot averaging, which is recommended for interpolation (the case N, = N ) [6]. In any case the N t-values
where x stands for some positive value and • for zero. The (N, X N , ) coefficient matrix for the normal equations to solve (1) turns out to be positive semi-definite, symmetric and to have a band width less than k , [6]. This makes the system suited for solution by Choleski factorization using De Boor's routines L2APPR, BCHFAC and BCHSLV. The same procedure is followed to obtain y ( u ) and z(u) to approximate the D,. and D. within tolerance ~ ~ and ½T: respective(y. In general this would lead to different numbers of con]
•
~;
6
State-of-the-Art
Computers in Indust~
ill
i
o
i
l
1
I ill 2
3
1tll I
J~J~____J_AL
0
2
!
3
4
S
6
7
: ...... :7
/ 3
fD .
__
!
[]
t
,
'il
/ /
[]
data points
•
curve points
L-]i'
j,
\
/
T× = lo i
L. . . . .
I ......
0
LENGHT
50
~,
=
10
d
100
UNITS
jJ :),
jl /~ //
k
¢
' i
~[
:i ~ /
/
L
i'
i
e
k
i
d
L
Fig. 1. Refitting a B-spline curve on a new uniform knot vector with a larger number of knots, while avoiding to get closer to the data points than necessary. Example projected onto the XY-plane. (a) 15 data points, B-spline curve, its control polygon and its knot vector to approximate the data within tolerance 7~, = 10 and T~. = 10 length units, with minimal number of control points (N~ = 6). Suppose that fits for the other curves imply that the number of control points must be increased to 10; then 10 abscissae are determined in the parameter range as described in the text. These and the corresponding curve points have been calculated and are shown in the figure. (b) B-spline curve with its control polygon that interpolates the 10 points evaluated on the curve of Fig. l(a). In the new knot vector the 10 abscissae used for interpolation are indicated. (c) Comparison between the curve of Fig. l(a) (indicated by arrows) and the refit of Fig. l(b). (d) Comparison between the curve of Fig. l(a) (indicated by arrows) and the B-spline curve that would be found by refitting the original data with 10 control points.
trol points for the x-, y- and z - c o m p o n e n t s , which w o u l d not produce a curve r(u) with three-dimensional control points. Therefore, in the i m p l e m e n tation the tolerance checks are carried out for all
three c o m p o n e n t s , for each N . . If o n e or m o r e of the c o m p o n e n t s deviate too m u c h from the data, N. is incremented for all three fits with x ( u ) , y(u) and z(u).
Cornputers in lnd~t~
J. Vergeest / Surface Fitting
At this point we have obtained a B-spline curve of order k,: U, = E i=1
that approximates ones row of data points in the least-squares sense, and in addition does not deviate from the data by more than half of the specified tolerance T.
3.3. Re-Representation of a B-Spline Curve on New Uniform Knots For each of the M rows of data points a B-spline curve ~)(u) can be determined in the way described above: N,,j
rj(u) = Z P, j B [ " ( u ) ,
j = 1 . . . . . M.
Note that the points P,j will in general not form a rectangular array. To stay within tolerance ±2 7", a different number N,, of control points may be needed for the different rows, depending e.g. on the arbitrary number of points contained in the rows. On the other hand, in the following stage of tensor product surface fit, we will need curves rj(u) on identical knot vectors and hence with identical N,; The obvious choice is to take N,m., = maximum of N~,, .... N,M, which implies that each N,, must be increased by (N, .... - N , , ) . There are several ways to do this. One way is to replace N,, by N,,,,. in the curve fit for r(u). However, if N,m.. > Nt for some j ( N / i s the number of data points in row j), the systems gets underdetermined. Even if N, .... ~< Nj, the newly fitted curve r(u) would have the disadvantage of being less smooth than necessary for the required tolerance (see Fig. l(c)). Therefore, we have chosen to keep the geometry of the curves r~(u) principally unchanged while increasing N , . This is accomplished in the following way. N,,,.... locations e,j on curve ~(u) are evaluated at parameter values t,, i = 1. . . . . N,m., e,, = ~/(t,),
where t, is determined in the parameter range of r(u) as described above, except that N (the number of abscissae) is replaced by N.m.. Then the locations e,j are regarded as N~..... data points to
7
which a B-spline curve of order k , with N u...... control points has to be fitted (which is in this case equivalent to interpolation). This too is performed in the way outlined above: the appropriate uniform knot vector is set up and Num~. abscissae are placed in the parameter range. The fitted curve will exactly interpolate the N, ..... data points and it can be verified that under the actual conditions of knots and t-values, the deviation of the curve from the point set r~(u) is small compared to the tolerance T, except for some extreme cases with large C k,,- ~ discontinuities. A detailed analysis would be required on to how ~(u) can be represented optimally under the given constraints. Note that increasing the number of control points from N~j to N. cannot be simply performed by knot insertion. This would in general produce a non-uniform knot vector, which is not allowed by the embedding SIPSURF system. Our re-representation is on basis of a new, but still uniform knot vector, for which in general an exact solution does not exist. Any underdetermined system (curve with N/< k . ) possibly encountered earlier in the scheme, will be solved here. If ~ < k . then the exactly interpolating B-spline curve ~(u) of order Nj is determined, on which N...~ locations are evaluated in the way just described, to which ~(u) of order k. with N.m.. control points is fitted.
3.4. Obtaining a Surface from the Fitted Curves We have now obtained M B-spline curves r~(u) of order k. with N.m., control points E,I each ( i--- 1 . . . . . Nu.,~ ; j = 1 . . . . . M ). Next we regard these control points as N,....... rows each consisting of M points. In principle the procedures described in the previous two sections are repeated, but now in the inverse direction. For every row i we determine by least-squares approximation the Bspline curve r,(v) of order k, that approaches the M points E,/ ( j = 1. M ) within tolerance ½T with the minimal number of control points N,,. By using exactly the same knot and abscissae conventions as for the curve fits described above and by reparametrizations, we will obtain N, ..... B-spline curves ~(v) with identical numbers (N,.m.) of control points P,j (i = 1 . . . . . N ...... ; j = 1 . . . . . N,,m.,). The number N , , will depend on the data, but it will be less than or equal to M and is the largest N,, encountered. Due to properties of
8
State-of-the-Art
Uomputers in Industry
Fig. 2. Surface fits with tolerance T = (0, 0, 0). (a) Data points and B-spline surface fit with degree = 1 (k,, = k,. = 2); (b) degree = 2; (c) degree = 3; (d) degree ~ 9; (e) overshoot in the fitted surface with degree = 9.
Computers in Industry
J. Vergeest / Surface Fitting
Table 1 Fits of a bicubic B-spline surface ( k , = k,, = 4) to three different p o i n t sets at different tolerances N u m b e r of d a t a rows
N u m b e r of data points
Tolerance
N u m b e r of coefficients
6 6 5 5 5 35 35 35 35
21 21 160 160 160 1775 1775 1775 1775
1.0 10.0 1.0 10.0 20.0 1.0 10.0 20.0 100.0
5X 6 = 5× 6 = 36x 5 = 28 X 5 = 20x 4 = 61 × 35 = 43 × 34 = 32 × 34 = 4×20 =
30 30 180 140 80 2135 1462 1088 80
Elapsed time a (seconds) SFIT01
(general)
1 1 10 7 4 293 150 87 12
(15) (15) (8870) (1200) (350) (-) (-) (23600) (955)
Elapsed times are as m e a s u r e d for runs on a Digital VAX-Station 3200.
tensor product surfaces, the control points P,/ are the (N,m,~ × N,,m,~) control points for the B-spline surface S(u, v) of order ( k , , k,,) defined on the former knot vectors for the u- and v-parameter that approximates the points E u in the leastsquares sense. As the curves ~)(u) deviate at most approximately ½T from the original data points D,j, while the curves r,(v) are closer than ~ T to the control points of the curves ~)(u), the surface S(u, v) approximates the original data within T, except for perhaps some extreme cases. In the current implementation there is no explicit distance check between D,j and S.
4. Implementation The fitting system has been implemented in standard Fortran as a computer program named SFIT0I. Interfacing with the interactive SIPSURF system proceeds, momentarily, through data base level operations. The user produces, somehow, a list of coordinate values of the data points in free format. The data may originate from a 2-D or 3-D digitizer, or be entered directly into a text file. The desired surface degree in u- and v-direction and the x-, y- and z-tolerance must be specified. Following the method described above, SFIT01 adds the new B-spline surface to the SIPSURF database. The user gets a report of which data points deviate maximally from the surface, and in which direction. It is the intention to arrive at a full integration of SFIT01 with the workstation graphics in the future.
4.1. Processing times" In the current implementation SFIT0I spends most of its time on finding the smallest knot sets, rather than on calculating the final control points. The program can be significantly optimized in the knot searching part. Still, as long as the number of data points is small, the total time to find the surface that passes the tolerance tests is short compared to the time needed for raster image generation on the particular work station. Table 1 lists the elapsed time in seconds to produce a SIPSURF B-spline surface from the data points. The degree in both directions is set to 3 in all cases. The time needed by SFIT01 (which uses De Boor's L 2 A P P R ) has been compared to the time spent by a similar program in which no account is taken of the specific local support properties of B-spline curves, but where instead a general least-squares solver is applied (routine E04FDF from N A G Library [10]). The programs have been run on a Digital's VAX-Station 3200 including a graphical display of 1024 × 768 pixels with 8 bit planes. Generation of a high-lighted color picture in a 256 × 256 pixel window typically takes 60 seconds. In the next section the possibilities of integration of SFIT01 with SIPSURF to enable interactive fitting will be indicated. As can be seen in Table 1, the number of coefficients found for the B-spline surface can be larger than the total number of data points when the tolerance is small. This is due to the fact that a row consisting of few data points is eventually fitted with a curve having as many control points
10
State-of-the-Art
('omputers in Industry
Fig. 3. Fits with B-spline degree = 2 to the data points of Fig, 2 for different values of the tolerance T: (a) 100 cm: (b) 50 cm: (c) 30 cm; (d) 20 cm; (e) 10 cm; (f) 1 cm.
Computers in Industry
a
J. Vergeest / Surface Fitting
11
b
Fig. 4. Highlighted image (left) and NC milled part (right; scale 1 : 20) from the model of Fig. 2(c). The symmetric other half of the car has been added to the model by using SlPSURVoperations.
as the maximum needed for any other data row (Section 3.3).
4.2. Example As an example, a set of very few data points (a total of 75 divided over 20 rows) has been taken from a simple four-view drawing of a car model (advertisement leaflet from the manufacturer). In each row only 3 to 5 data points were determined on the left half of the car body. Figure 2 shows the data points and fits to the data at zero-tolerance with B-spline degree 1, 2, 3 and 9. For degree 1, 2 and 3, B-spline surfaces with 20 × 5 control points were obtained, for degree 9 one with 20 × 10 control points. The gross shape of the car body is reproduced with one single B-spline surface. The fit with degree 2 seems optimal, as it features the fewest overshoots (curling self-intersecting surface). To remove such overshoots, the user may add a few data points to guide the surface fit. Alternatively, the tolerance can be increased. In Fig. 3 the results of fits with degree 2 for different tolerance values are shown. B-spline surfaces of 20 × 5 control points are found for tolerances up to 20 cm (length of car is approximately 400 cm). For larger tolerances fewer control points are sufficient. Figure 4 shows an NC milled model out of foam material (after adding symmetrically the other half of the car body using SlPSUaF oper-
ations " G l u e " a n d / o r "Connec" for automatic G 2 blending). The question arises of what the user should do to further modify and detail the shape. A standard option of SlPSUaF is control point editing (designation with a mouse and typing delta coordinate values). However, the control points found by fits are not easily interpretable generally (Fig. 5). It will be investigated in the future what the possibilities are of moving and adding data points (rather than control points) in an interactive way.
5. Conclusions and Further Research
5.1. Conclusions A method has been presented to generate a single B-spline surface on the basis of (possibly very few) data points. Smoothness and closeness to the data is controlled by simple tolerance values. The program SliT01 automatically finds the B-spline surface with the smallest number of control points to meet the tolerance requirement. We have illustrated that it would not have been possible to find those control points in a classical trial-and-error procedure, directly by the user. This implies that further control point editing after the fit is not practical in general. On the other hand, since the fitting procedure is sufficiently fast for
12
State-of-the-Art
Computers in lndustrv
-: j )
Fig. 5. The control points found by SF~T01are often not suited for direct further modeling by the user (shown are the control points of the surface of Fig. 2(c)).
less than about 300 data points, interactive modifications, addition and deletion of single Or groups of data points seem feasible in implementations on present workstations. 5.2. Further Points o f Interest
There are alternatives to some of the choices made for the fitting procedure. Moreover, several improvements and extensions to the currently implemented system are possible. Some of the issues are: (1) Most of the computation time in SFIT01 is spent on finding the simplest possible B-spline surface, not on fitting its control points. Much time could be safed if, instead of simply incrementing N, by one (Section 3.2) a more direct search strategy (e.g. a binary search) would be followed to obtain the optimal N, for the curve fits. (2) The maximal deviation of the surface from the data points is the criterion for SFIT01 to accept the fit or to increase the number of B-spline control points. A direct way to check against this criterion is possibly with surfaces that minimize the L ~ norm (largest deviation). Such fits can be implemented with N A G Library routine E02GCF, but these are as slow as the general solver (E04FDF, see Section 4.1). We could not find an L ~ minimizer for B-spline curves that could approach the speed of the De Boor's L 2 (least-squares) algorithm, as implemented in SFIT01. (3) SFIT01 can be extended to closed and doubly closed B-spline surfaces (the three other topologies of SIPSURF, see Section 2.1). This would be
very appropriate if the data points are known to be located on closed contours. For closed (periodic) curves, the matrix of size iN,, × N ) in Section 3.2 will contain a few more non-zero elements in the top-right a n d / o r in the left-bottom corners, causing the loss of the band structure of the normal equations matrix. The system can still be solved with the general N A G Library routine E04FDF. To take advantage of the (still highly) structured matrix, algorithms due to Dierckx can be applied [7]. (4) The parametrization of any fitted curve is independent from the geometry of nearby other curves (except the size of the knot vector). Suppose that the M curves found in Section 3.4 happen to be well alligned or even parallel in a plane. If one of those curves was obtained from a row of data points with a large concentration of points at a particular place, then there will be also a concentration of control points of the fitted B-spline curve near that place. If nearby data rows do not exhibit such a concentration, the surface as a whole, found in Section 3.4 may be distorted. For some applications this may just be the desired effect. In many cases, the surface should depend on the geometry of the intermediate curves, not on their parametrizations. Parametrization independence can be obtained by evaluating points on a curve at intervals of equal curve length and then refitting the curve to those points. This could be an alternative for the re-representation described in Section 3.3. For fits to data on closed contours, the reparametrization should even allow shifts of
Fig. 6. Unwanted twists and folds in the surface if the parametrization does not account for irregularly spaced data points.
J. Vergeest / Surface Fitting
Computers in Industry
the begin/end locations of the curves to avoid unwanted twists in the resulting surface. Figure 6 gives an example of a surface obtained by SFIT0I from uncorrected photogrammetric contour data. (5) The curve fits in SFIT01 are unweighted. If there is knowledge about relative weights of the data points, it is straightforward to set of w s in (1) accordingly. It is not clear what the effect of the weights should be in the further stages of the procedure. One possibility is to propagate the data point weights to the relevant curve control points for the transverse fits. (6) The modeling capacity of SFIT0I when integrated with the SIPSURr surface modeler needs further investigation. When raster image display is sufficiently fast, will then the user be able to obtain shapes with all necessary features? Unwanted fluctuations should be damped by adding or removing data points, or by changing their weights. (7) SFIT0I can be significantly speeded up for operations that effect only one data point or one row of data points (these are the operations most likely to occur interactively). The curve fits on the other rows then need not be redone in many cases. The transverse fits must all be repeated, but a large part of the matrix equations remains unchanged. (8) User interface and display graphics should be optimized for interactive shape design. Matters of interest are movement of data points in 3-D, fast coarse image generation with automatic refinement until the next user interrupt, storage of a number of recent images including the data from which the user operations causing those images can be deduced.
Acknowledgements The author wishes to thank Wim Sybranda for making available preliminary design data, Rinus
13
Grabijn for preparing illustrations and Bram de Smit for the production of the NC-model.
References [1] E. Andersson, R. Andersson, M. Boman, T. Elmroth, B. Dahlberg and B. Johansson, "Automatic construction of surfaces with prescribed shape", Comput. Aided Des., Vol. 20, No. 6, 1988. [2] R.E. Barnhilt, "Representation and approximation of surfaces", in: J.R. Rice, ed., Mathematical Software 111, Academic Press, New York, 1977, pp. 69-120.' [3] R.E. Barnhill, "'A survey of the representation and design of surfaces", 1EEE Comput. Graphics AppL, Vol. 3, No. 10, October 1983, pp. 9-16. [4] W. Boehm, G. Farin and J. Kahmann, "A survey of curve and surface methods in CAGD", Comput. Aided Geom. Des., Vol. 1, No. 1, 1984, pp. 1-60. [5] J.J. Brock, "CADAMP, Computer Aided Design and Model Production", Proc. DECUS '84, Digital Equipment Corporation, Marlboro, MA, 1984. [6] C. de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978. [7] P. Dierckx, FITPACK User Guide, Part 1: Curve Fitting Routines, Report TW 89, Dept. of Computer Science, Katholieke Universiteit Leuven, Belgium, 1987. [8] J.G. Hayes, "Fitting surfaces to data", in: R.R. Martin, ed., The Mathematics of Surfaces 11, Clarendon Press, Oxford, 1987, pp. 17-38. [9] N.J. Lott and D.1. Pullin, "Method for fairing B-spline surfaces", Comput. Aided Des., Vol. 20, Nr. 10, 1988, pp. 597-604. [10] NAG, The NAG Fortran Library Manual--Mark 12. Numerical Algorithms Group, Oxford, 1987. [11] J.S.M. Vergeest, "Connecting arbitrary surfaces under geometric constraints", Computers in lndustlT, Vol. 8, No. 1, 1987, pp. 3-12.