Smooth-surface

approximation and reverse

engineering
B Sarkar and C-H Menq

representation of a scanned object with multiple

The paper identifies the steps involved in the reverse- surfaces. When a surface containing multiple features
engineering process. The procedure begins with the (cylinders, cones, spheres, sculptured surfaces etc.)is
division of the whole array of measurement data points scanned, it is imperative to detect the shape changes
into regions, according to shape-change detection. In along the boundaries of these features. Scanned points
each region, points are parameterized, and knots are are divided between different regions on the basis of
selected. Smooth parametric surface approximation is the location of the boundaries, and a surface can be
obtained by the least-square fitting of B-splines. Nonlinear fitted to the points belonging to each of these regions.
least-square minimization is applied for parameter B-spline surface representation is used in this paper. It
optimization with simple bounds on the parameter has many advantages4:
values. The objective function minimized is the explicit
error expression for the sum of the squares of error • One surface is considered as a single element, which
values at the data points. avoids inner constraints, such as the specification of
derivatives and twists.
reverse engineering, CAD~CAM,surfaceapproximation • B-spline basis functions have local support, so that
modifications affect a limited area.
• A relatively small number of points is required for
complicated geometries.
With the advent of high-resolution laser scanning
machines and computer-controlled measurement In each region, an explicit expression for the sum of
machines, designers can now create and/or modify the the squares of the error values at the data points is
design of an object based on a manufactured part that obtained. This expression entails the evaluation of
might have been modified from the existing design B-spline basis functions at the parameter values. Simple
while in production. The industrial jargon for this bounds are applied to the parameter values of the
process is reverse engineering. In the die and mould points that are internal to a region. The parameter
industry, the existing design is often modified on the values at the points on the boundary are fixed. The
shop floor, depending on the production requirements. error expression is then minimized. Once the parameters
The product geometry changes with such modifications. are optimized, least-square fitting is applied. This
The idea behind reverse engineering is to retrieve method precludes repetitive fitting, and, although the
this new geometry from the manufactured part computation time can be long, depending on the
by scanning, and to modify the existing CAD number of points, it leads to a full surface approximation.
database. The coordinates and the normal vectors of The designed surface is then subdivided into a
the part surface, with respect to a given datum, are user-defined number of patches to achieve a prescribed
the entries in a CAD database that are often modified tolerance level. Figure 1 gives an overview of the whole
on the basis of the output of the reverse-engineering procedure.
techniques. These techniques mainly involve parametric
surface approximation.
A few papers have been published on the problem
of approximating 3D data by parametric surfaces 1-3. STEPS IN REVERSE E N G I N E E R I N G
However, these papers describe the approximation In this section, the steps involved in the generation of
procedure for a single isolated surface, and the surfaces are described in detail. The input data consists
optimization of parameters involves repetitive fitting. of points in R 3, and the generated surface may not
The objective of this paper is to provide smooth deviate from these points by more than a given
tolerance. These points may be measurement points
Department of Mechanical Engineering, Ohio State University, obtained from a coordinate-measurement machine
Columbus, OH 43210, USA (CMM) or a laser scanning machine, or they may be
Paper received: 27 February 1990. Revised: 15 November 1990 obtained from the existing CAD database.

volume 23 number 9 november 1 9 9 1 0010-4485/91/090623-06 © 1991 Butterworth-Heinemann Ltd 623

 Scansurface ~ - ~ Detectboundary ~ - ~ Divideintoregions I

Objectmodel

Surfaceintca~scction

Figure 7. Flowchart of reverse engineering Figure 2. 3D plot ot scanned points on turbine-blade die
model
,,'1
Boundary detection
Boundary detection is the primary step in reverse
engineering wherein the total array of input data is
divided into several regions. Each of these regions is
devoid of any sharp change in shape, and consists of
a single feature. This procedure is sometimes automatically
performed in laser scanning machines s. In the case
where shape detection is not inbuilt into the software
of the measurement machine, one can apply the edge
operators used in computer-vision techniques• A
difference formulation of the edge operators, such as
the Laplacian Gaussian operator V2G, can be applied
to detect a change in shape 6. Here, the z coordinate
at each of the points scanned takes the place of
grey-level intensity at the pixels of an image. Edges are
identified by detecting the zero crossings, and nodes

Ga gR1111
on neighboring edges are connected to obtain the
boundaries. The efficiency of this algorithm depends
on the particular edge operator selected and the nature
of the data.
The following example provides an insight into the
boundary-detection method. A turbine-blade die wax
model was scanned for 3000 (60 x 50) points on a
Sheffield Cordax RS-30 CMM. The edges were detected
by convoluting the z coordinates of the measurement
points with the V2G operator using a 6-point mask.
Figure 2 shows the primal sketch of the scanned points,
and Figure 3 shows the3D plot of the points with zero
crossings at the boundaries after the use of the V2G
operator. The blade model was then divided into
regions, as shown in Figure 4. Points belonging to each
region were identified• Region 1 was the blade profile
surface, and a B-spline surface was fitted to it.
Region4
Figure 4. Detected regions in turbine-blade die
Parameterization
Parameterization is the process of assigning parameter Along each track, chord-length parameterization is
values to the measurement data points in each region. applied to assign u values. The same procedure is
Given an ordered set of points Pij(xii, Yii, zij), the repeated for tracks defined by u = 0 and 1 to assign
parameter value (ui, vi) can be assigned by chord-length the v values. The v parameters for the interior points
parameterization 1'7 or equal-increment parameter- are assigned by linear interpolation•
ization 8. Chord-length parameterization is used in this
paper, because it bears a close relationship to the actual
Selection of knots
spacing of the data points. In the scanning of a surface,
points are generally taken along lines aligned with one The designer chooses the number of knots, or, in other
of the machine axes. These lines are called 'tracks'. words, the number of spline segments in both

624 computer-aided design

parametric directions, and the parametric values at and ~.i. Mi(u) is 0 everywhere except in the range
these knots. The number of segments clearly depends 2i_" < u < 2i. Similarly, let the B-spline basis function
on the number of data points and the prescribed Nj(v) be defined in the v direction; it is nonzero in
tolerance value. As the number of knots increases, the 6i_ 4 < v < 6i. The composite basis function for surface
sum of the squares of the individual errors decreases. fitting is a product of the B-splines Mi(u) based on the
The parametric value at these knots, on the other knots ]i and the B-splines Nj(v) based on the knots 6i,
hand, is governed by the geometric factors, such as as shown in Figure 5. The general bicubic spline surface
the position of the inflection points or the nature of is represented in the form
variation of the curvature. The designer should consider h+4k+4
having more knots at those portions of the curve or s(u, v) = ~ ~ ciiMi(u)Ni(v) (5)
surface where rapid variation of the curvature is implied i=1 j=1
by the data. Although adaptive generation of knots The observation equations for the least-squares method
and automatic parametric-value assignments can become
reduce the intuitive work required in this step, the h+4k+3
designer can, nevertheless, obtain some flexibility when ~, ~ c~iMi(u~)Nj(v~) = Xr (6a)
selecting the knots interactively. i=1 j = l
h+4k+4
Least-square fitting ~, c~Mi(u~)Nj(v~) = y~ (6b)
i=1 /=1
: Least-square fitting follows the parameter optimization h+4k+4
procedure in the sequence of steps for reverse E E c~jMi(ur)Nj (vr) = Zr (6c)
engineering. However, as the formulation of the error i=1 j = l
expression for optimization requires the mathematical where r = 1, 2,..., ran.
background of least-square fitting, this step is described Equation 6a can be written in matrix notation as
first.
Hayes eta/. 9 describe the mathematical background TC ~ = X (7)
for cubic B-spline fitting to scattered data. The objective The solution of this equation by means of the
of this step is to find the functions x = flu, v), pseudoinverse is
y = g(u, v), z = h(u, v) for a surface. Mathematically,
in least-square surface fitting, one aims at approximating C x = (TTT)-ITTX (8a)
a set of data points (xq, Y~i,z~j), where i = 1, 2,..., m Similarly, the solutions for Equations 6b and 6c are
and j = 1, 2 , . . . , n (ran being the total number of
measurement points) by cubic B-splines with selected C Y = (TTT)-ITTy (8b)
knots. The rectangle ranging over [0,1] in the u C z = (TTT)-ITTZ (8c)
direction and [0, 1] in the v direction contains all the
.data points of a particular region. The rectangle is Here, T is a matrix with m rows and (h + 4)(1 -t-4)
subdivided into smaller regions by selected knots that columns, and C x, C v, and C z are vectors with
satisfy (h + 4)(1 + 4) elements. X, Y and Z are vectors with
0 = ~'0 < ~1 < ~2 < "'" < ~'h+l = 1 (la)
0=60<61 <62<...<6k+1=1 (lb)
where h and k are the number of knots in the u and Mi (u) Ni (v) surface
v directions, respectively. The knot sets are augmented
by 16 otherwise arbitrary additional knots satisfying III I I% I I | I I II
~'--3 < ~'--2 < ~'--1 < ~0 = 0 (2a)
1 = 2h+1 < 2h+2 < 2h+3 < 2h+4 (2b)
6_ 3 < 6 _ 2 < 6 _ 1 < 6 0 = 0 (3a) ~j--1 H - - H
~j--2
1 = 6k+ 1 < 6k+ 2 < 6k+3 < 6k+4 (3b)
I ~j--3
The 4th-order (cubic) B-splines can be evaluated by
~i-4
the following recurrence relationship:
v
(U -- ,~i_4)M3,i_ l(U) -{- (~i - LI)M3,i(U)
M4,i(u) =
2 i - - 2i__ 4
(4a)

1
MI'i(u) -- /]'i- ~i-1 '~i-1 ~ u ~ '~'i
(4b) ~--3 ~-2 ~-1 0 ~i-4 ~i-3 ~i--2 ~i-1 ~i
= 0 otherwise U .~
Let the B-spline basis function Mi(u) be defined in the Figure 5. Rectangular regions created by knots
u parametric direction with knots 2i_4, ~i-3, 2i-2, 2i-1 [The basis function is nonzero over 16 adjacent regions.]

volume 23 number 9 november 1991 625

m n elements. The element of the matrix 1" in the rth The expression for the total error value becomes
row and [(j - 1)(h + 4) + i]th column is Mi(u~)N,(vr). mn mn
The elements of the vectors C x, C Y and C z ' are E(ur, vr)r=l,2 ........ = ~, ( e y = ~, (e~)2 + (eY)2 + (e~)2
arranged in the sequence i = 1, 2.... , h + 4 for fixed j, r--1 r=l

wherej = 1, 2,..., k + 4. Once the optimal parameters = XTPX + YtPY + ZrPZ (11)
are determined, the problem is reduced to that of
finding the coefficients c~, c~, c~ in the observation The average error is defined by
Equations 6a-c. The Householder transformation of [m~ ]1/2
matrix 1" can handle the solution better in the ear = ~, ( e r ) 2 / m n (12)
r=q
rank-deficient cases and when 1" is very ill conditioned,
but the computation time is doubled. F(ur, vr) is the objective function for the nonlinear
least-square minimization problem with u~, v r,
r = 1, 2 . . . . , mn, as the variables. Simple bounds are
Parameter optimization
imposed on the parameter values at the points internal
Hoschek 2 describes a linearized approach to the to a region, and constraints are imposed at the
optimization of the parameters with the object of boundary points of the region of interest. That is, u i = 0
making the error vectors normal to the fitted surface. or 1 and vi - 0 or 1 at the edges of the parametric
In his method, each error vector is projected onto the space, and 0 < u i < 1 and 0 < vi < 1 for the internal
tangent at the point on the fitted surface that points. This nonlinear least-square problem may be
corresponds to the measured point. The projected solved by the modified Levenberg-Marquardt method
length after normalization yields the change in and an active set strategy with numerically computed
parameter values. Pratt I mentions a similar algorithm, Jacobians. In all the cases tried so far, the average error
but he solves nonlinear equations in making the error eav decreases significantly after only a few iterations.
vectors normal at the respective points. In both of The following example illustrates the surface-fitting
these methods, the parameter values are sought after procedure, and compares the above-mentioned
a specified number of iterations, and the surface is optimization method with Hoschek's method with
refitted during every iteration. respect to computational efficiency. The data points
This paper uses a different method for parameter for this example were obtained by scanning a wax
optimization. The difference between the measured model of a surface patch designed at CATIA (Dassault
data points and the corresponding points on the fitted systems) and machined on a Cincinnati Milacron
surface is minimized in the least-squares sense, and no machining center. The scanning was performed on the
refitting of the surface (which is computationally costly) Cordax CMM. The data consisted of 100 points: ten
is required. An explicit expression for the sum of the points on each of the ten section curves. The original
squares of the error vectors is obtained, and it is then surface was designed with 4th- and 5th-order splines
minimized, with the parameters as the variables, by in the two parametric directions.
the use of a nonlinear least-squares minimization Table 1 shows the total error value before and after
technique. optimization for different numbers of knots selected,
The error components in x at each point are obtained with the CPU time used. The optimization procedure
in matrix form: is carried out for up to ten iterations. It may be observed
that there are several ways to design a surface for a
Ex = x -- T ( T t T ) - I T t X = [I -- T(TTT)-ITT]X = PX prescribed tolerance. For example, a tolerance value
(9a) of, say, 5.0 x 10 -3 mm can be achieved in various
ways. Starting with a higher number of knots may
where P -- [I -- T(TTT)-ITT]. The elements of vector require fewer iterations (four iterations for I x 1 knots,
Ex are the x components of the error at every point three for2 x 2 a n d 2 x 3, and one for3 x 4 a n d 4 x 5
e;, r = 1, 2 , . . . , m n . cases). However, the CPU time increases with the
The y and z components of the error are similarly increase in the number of knots for a given number
obtained in the matrix form as of iterations. On the other hand, if one starts with fewer
knots, more iterations are required to attain the same
Fy = py (9b)

Ez = PZ (9c) Table 1. Variation of error with number of knots and CPU

mn
usage
(eX)2 ---- EtEx = XTPtPX = XTPX (10a)
r=l Number of Number of Average Average CPU time
segments in segments in error before error after
as p t = p and p2 = p. Similarly, u direction v direction optimization optimization
mm mm s x 10 ~
n'In
(ep2 = yTpy (10b) 1 1 1.412 0.111 1.79
r=l 2 2 0.668 0.031 2.56
2 3 0.647 0.029 3.08
mn 3 4 0.208 0.014 4.49
(el) 2 = zTpz (10c) 4 5 0.169 0.005 6.96
r=l

626 computer-aided design

error tolerance. The computation time increases in this
case, also. Therefore, given an error tolerance value,
the optimum number of knots needs to be determined.
This problem has been studied at length by Sarkar and
Menq 12.
So that the suggested optimization method could
be assessed, Hoschek's method was implemented using '\
B-spline basis functions. The method in this paper was
found to be computationally efficient in all the cases
studied so far. In the above case of 3 x 4 knots, after
4493 s of CPU time, the method achieved an average
error value em of 0.014 mm after ten iterations starting
with an error value of 0.208 mm. Hoschek's method
produced an error value of 0.0195mm after 815
iterations in the same amount of CPU time 12. Although
the CPU time required per iteration for Hoschek's
method is much smaller than that for the method in
Figure 7. B-spline surface fitted to measurement points
in Region 1 of Figure 4
this paper, the number of iterations required to achieve
a particular error tolerance nevertheless more than
compensates for the total CPU usage. Figure 6 shows
the surface fitted to the measurement points with three gone through. In this case, the number of patches, or,
knots in the u direction and four knots in the v direction. in other words, the number of segments in the u and
Finally in this section, the results of the optimization v parametric directions, would have to be increased,
algorithm applied to Region 1 of Figure 4 are described. leading to the introduction of new knots and
The points in this region were assigned parameters by subdivisions. Clearly, the number of subdivisions
the use of chord length. The turbine-blade model depends on the error tolerance value. The main idea
consisted of 30 x 11 patches in the original CATIAmodel, is that the design starts with a rough form, and it
on the basis of which it was NC machined. This region converges to the final form after undergoing subdivisions.
was modelled from 120 (15 x 8) measurement points Therefore, starting with the selection of knots, achieving
with 8 x 4 knots in the u and v directions, respectively. the right tolerance value is an iterative procedure.
The results of the optimization after five iterations were:
average individual error before optimization Surface intersection
= 0.004 01 mm
In this step, the boundary between two surfaces fitted
average individual error after optimization
in two adjacent regions can be recomputed using a
= 0.001 83 mm
surface-intersection algorithm. Accurate calculation of
CPU time taken = 16 320 s
the boundary may be essential if the surfaces modelled
Figure 7 shows the blade surface reconstructed using by reverse-engineering techniques are to be NC
B-splines as a dense net of orthogonal parametric lines. machined. A number of algorithms are available for
the parametric surface-intersection problem. An
Surface subdivision overview of surface-surface intersection methods has
been presented by Pratt and Geisow 1°. The method
It is possible that the prescribed tolerance value could used for this work falls into the category of lattice
be still not achieved after the above steps had been evaluation. Intersection points between parametric
curves on the primary surface and the second surface
are determined. After these points have been obtained,
the usual problem is to sort them out in a presumable
order so as to define the whole intersection curve. In
the problem in this paper, the surfaces that intersect
lie topologically adjacent to each other, and can thus
intersect only once near the predetermined boundary
between them. This fact also helps in the incorporation
of the initial guesses for the numerical solution of the
problem. Once determined, the intersection points are
connected with a B-spline curve 11,
If p(u, v) and q(s, t) are two adjacent parametric
surfaces, the distance between two arbitrary points
lying on the two respective surfaces is
d(u, v, s, t) = r(u, v) - q(s, t) (12)
If intersection points between the u parametric curves
of surface p(u, v) and surface q(s, t) are sought, the
Figure 6. B-spline surface fitted to "100 data points problem is reduced to three variables. The distance

volume 23 number 9 november 1991 627

vector becomes Ohio State University, USA, for supporting this research
project.
d(v, s, t) = p(ul, v) - q(s, t) (13)
The distance vector d(v,s,t) is minimized in a
REFERENCES
least-squares sense by the use of the nonlinear
Levenberg-Marquardt minimization algorithm. Fast 1 Pratt,M J 'Smooth parametric surface approximations
convergence is ensured by the fact that there can be to discrete data' Comput. Aided Geom. Des. Vol 2
a single intersection point between the parametric (1985) pp 165-171
curves of one region and the surface in the next 2 Hosehek, J 'Intrinsic parametrization for approxi-
adjacent region, and the parameter values are close mation' Comput. Aided Geom. Des. Vol 5 (1988)
to 1 or 0, depending on the assignment of the pp 27-31
parametric directions in these regions.
3 Lyche, T and Morken, K 'Knot removal for
parametric B-spline curves and surfaces' Comput.
Aided Geom. Des. Vol 4 No 3 (1987) pp 217-230
4 Riesenfield, R 'Application of B-spline approximation
CONCLUSIONS
to geometric problems of computer aided design'
The steps for surface approximation in reverse PhD Thesis Syracuse University, USA (1973)
engineering are explicitly described in this paper. A
method has been developed to divide the array of 5 Miyoshi, T 'Mold production/machining system
points obtained from scanning an object with multiple based on non-contact digitizing system' OI(ADA
features into several regions. Each region is fitted with Vigitizer Technical Information
a surface. B-spline representation of curves and surfaces 6 Marr, D Vision Freeman, USA (1982)
has been found to be extremely advantageous. Simply
by increasing the number of segments, one can achieve 7 Hadley, P J and Judd, C J 'Parametrization and
a desirable error tolerance. The optimization procedure shape of B-spline curves for CAD' Comput.-Aided
used in this paper is also novel, in the sense that an Des. Vol 12 No 5 (1980) pp 235-238
explicit expression for the error value is obtained, and 8 Faux, I D and Pratt, M J Computational Geometry
it is then minimized by a nonlinear least-squares for Design and Manufacture Ellis Horwood, UK
optimization technique with the assigned parameters (1980)
as variables. Only a few iterations can yield a very
9 Hayes, l G and Halliday, l 'The least-square fitting
prudent reduction in the error. The method has been
of cubic spline surfaces to general data sets' J. Inst.
found to be computationally efficient. Last, the fitting
Math. Applic. Vol 14 (1974) pp 89-103
procedure has to be applied only once after the
optimization of the parameters, as opposed to the 10 Pratt, M J and Geisow, A D 'Surface/surface
repetitive fitting used to date. The authors feel that the intersection problems' in Gregory, J A (Ed.) The
sequence should be quite useful for implementation in Mathematics of Surfaces Clarendon Press, UK (1986)
a computer, regardless of the method used for pp 117-142
representing surfaces.
11 Golden, M E 'Geometric structural modeling: a
promising basis for finite element analysis' Comput.
& Struct. Vol 10 (1979) pp 347-350
12 Sarkar, B and Menq, C H 'Parameter optimization
ACKNOWLEDGEMENT in approximating curves and surfaces to measure-
The authors would like to thank the US National Science ment data' Comput. Aided Geom. Des. (accepted
Foundation ERC for Net Shape Manufacturing at the for publication)

628 computer-aided design