- Email: [email protected]
A surface-lofting approach for smooth-surface reconstruction from 3D measurement data
Computers in Industry 34 (1997) 73-85
A surface-lofting
approach for smooth-surface reconstruction from 3D measurement data
Chih-Young Lin, Chung-Shan Liou, Jiing-Yih Lai Mechanical Engineering Department,
*
National Central Uniuersity, Chung-Li 320, Taiwan, ROC
Received 18 September 1995; revised 1 May 1996; accepted 27 August 1996
Abstract A reverse engineering method for sculptured surface modeling is proposed in this work. It includes a curve fitting algorithm and a surface lofting algorithm, where the former is used to fit each row of measurement data to obtain a B-spline curve and the latter ifs used to fit the section curves into a lofted surface. The above two fitting algorithms are essentially one-dimensional, whereas the commonly used B-spline surface fitting is two-dimensional. This makes the proposed method more effective in temis of computational efficiency and CPU memory required. Both computer simulation and experimental results indicate that I:he proposed method can efficiently be used for the reverse engineering of complex object shape. 0 1997 Elsevier Science B.V. Keywords: Reverse enguieering: Surface lofting; Curve fitting
1. Introduction Reverse engineering is a technology to establish CAD models from samples, prototypes, molds or manufactured parts by digitization. It can be used as a design tool for extracting the design concept of an existing product by analyzing its CAD model. The integration of rev’erse engineering and the CAM technology can even replace conventional copy milling method, commonly used to fabricate products from clays or wood models. There have been considerable studies with respect to fitting sculptured surfaces to known 3D measurement data [l-4]. Many problems still exist due to the complexity and variety of the shape of practical objects. Some of these are
* Corresponding author.
such as computational problem, the difficulty in processing enormous amount of data, the accuracy problem in created surface model, the difficulty for complex shape modeling, etc. The objective of this work is to present a surface lofting approach for the reverse engineering of complex shape. The problem of fitting parametric surfaces to known measurement data has been studied extensively. Numerous investigators have focused on the development of various surface optimization schemes for accurate modeling of sculptured surfaces. Some of the numerical methods proposed are such as the Newton-like approach [l], the first-order Taylor series expansion method [2] and the LevenbergMarquardt method [3]. These methods can generally yield results with comparable surface accuracy for an isolated sculptured surface. However, the Newton-
0166-3615/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved. PII SOl66-3615(96)00082-6
74
C.-Y. Lin et al. / Computers in Industry 34 (19971 73-85
like approach [l] and the first-order Taylor series expansion method [2] are computationally less efficient than the Levenberg-Marquardt method 131. But, the Levenberg-Marquardt method requires a huge CPU memory, which limits its application to small amount of measurement data. The main drawback of the above surface fitting methods are that when the shape of the object becomes complex, none of them can yield effective CAD models. The maximum dimension used in the Levenberg-Marquardt method is (mn) X (mm), where m and n are the number of measurement points along the x and y coordinate axes, respectively. Both the computer time and CPU memory would increase substantially when the number of measurement points increases. Sarkar et al. [4] extended the Levenberg-Marquardt method by introducing a boundary detection scheme to divide the compound surface into several regions, each of which was fitted by a surface patch individually. However, how adjacent patches were connected was not addressed. The robustness issue of the proposed boundary detection scheme was not investigated also. A computer-aided design technique called ‘skinning’ or ‘lofting’ could be employed for the reverse engineering applications. This technique interpolates a family of curves defined in one of the surface’s parametric directions to form relatively complex surfaces. Tiller [5] and Woodward [6] gave a detailed description of constructing a lofted surface with Bsplines as well as the necessary background about interpolation. Lin et al. [7] presented a B-spline Boolean-sum approach for surface interpolation. A bivariate B-spline surface was constructed from an adequate superposition of univariate B-spline interpolation by the Boolean sum operation. The procedure was the same as that of Coons-Gordon surfaces. A possible way of extending the above approach to NURBS is also discussed for the approximation of Coons-Gordon-type surfaces. Filip et al. [8] pointed out several problems with using tensor-product B-spline surfaces in lofting, such as difficulty in lofting non-planar curves, how to determine appropriate knot vectors and the ‘crease’ problem. The crease problem arises when the surface normal does not change continuously, which results in a surface that is not visually smooth. They proposed a knot vector function, instead of a tensor-
product form where a single u knot vector is specified for every u value, such that the problem of lofting non-planar curves can be resolved. A procedure was also described for reparametering a G’ curve into a C’ curve so that the crease problem can be eliminated. An important characteristic of a set of curves capable to be lofted is that they should be comparable in degrees, number of spans and knot vector form. Cohen [9] proposed an Oslo algorithm for adding new knots. Associated algorithm for attaining new control points was provided also. Boehm [lo] developed another algorithm for adding knots and calculating new control points. Such an algorithm is actually a special case of the Oslo algorithm, but is computationally more efficient. Cohen et al. [I I] proposed a degree-raising approach for adding flexibility to a set of B-spline curves.
2. Problem statement The primary limitation in the available surface optimization methods [l-4] is that they cannot resolve the problem of complex shape modeling. When the measurement data increases, the substantial increase in computational time and CPU memory makes these methods ineffective. Moreover, a single isolated surface is often not good enough to describe the shape of the object surface. It is generally necessary to divide the object surface into a number of small patches such that the surface optimization algorithm can be applied. The objective of this work is to develop a method incorporating a curve fitting algorithm and a surface lofting algorithm for the reverse engineering of complex physical parts. The curve fitting algorithm is employed first to fit each row of measurement data, which yields a set of parallel B-spline curves. A least-squares algorithm is employed to optimize the control points and parameter values for the B-spline curves, so as to minimize the average deviation between each row of measurement data and the curve fitted. The surface lofting algorithm is then employed for interpolating the above parallel curves to construct a lofted surface. The proposed surface-lofting algorithm is different from those described in the literature [5-81. The
C.-Y. Lin et al. /Computers in Industry 34 (1997) 73-85
main differences are explained below. Firstly, for the CAD application, fhe describing curves are designed in the CAD system and the lofted surface is expected to pass all describing curves. For the reverse engineering application, however, the created surface would distort if all section curves are to be passed. This is because the available surface-lofting algorithms are not suitable for curves which are densely distributed. A fittin,g algorithm is developed to fit the section curves into a smooth surface. Secondly, the surface accuracy is an important issue in surface reconstruction. The lofted surface must be able to describe all points on the part surface. A reparameterization process :is employed to match this goal. Finally, sharp edges and comers appear frequently in the reverse engineering problem, which may lead to unexpected oscillation on the created surface. A detail discussion o-F the above issues will be given and methods will be provided to solve these problems.
of U. The basis functions Ni,k(~) are related to the knot vector U, and are defined as [12]: 1 N,,I =
ifx;
i 0
N, (Uj =
fitting algorithm
The curve fitting problem has been extensively studied and various computational algorithms are available in the literature [1,2,4]. However, for the completeness of thi;s paper, the algorithm used in this work is briefly described. Each array of measurement points (along the x direction) is modeled by a B-spiine curve. The measurement points are densely distributed such that the cross sectional profile of the part surface can truly be described. An opeln cubic B-spline curve is used, which can handle b’oth smoothness and local control property for each array of measurement data to be fitted. The position vector P(U) of a B-spline curve along the parameter u can be expressed as P+l
2
1,
(1)
where B, are the position vectors of the control points, p + 1 is the number of the control segments, P&(U) are normalized B-spline basis functions and urnin and s,,, are ,the lower and the upper margins
(2)
otherwise, (“-xi>N,,k-l(u)
r,k
x;+!c-1
+
(Xifk
--xi
-
wi+wG4
xi+
k -xi+
3
(3)
1
where xi are the knot values of the knot vector X. The normalized form of a kth order open uniform knot vector X is given as [12]: X =
[ x1 x2
. . . xkXk+l
xi = 0, 1 I ilk;
... x,+1 x,+2 ... x,+,‘+l
xi =
k+l
3. Curve
15
1f (4)
i-k p+2-k’
xi=1,p+2~i
The curve fitting algorithm is employed to optimize the position of the control points Bi and the parameter value u, as given in Eq. (11, corresponding to each measurement point. The other parameters, e.g., knot vectors and the number of control segments, in the B-spline model are determined beforehand. The average error eav is defined as the average of the sum of the squares of the fitting errors at each data point, i.e. e a” ZZ
;)!,
Cej)’
’
(5)
Je,=Iq-P(uj)I,
(6)
where m is the number of the measurement points, ej represents the shortest distance between the jth measurement point 5 and the curve fitted. The parameter optimization problem is formulated as the minimization of the error function eaV which yields an optimum set of solutions for Bi and uj. An initial estimation of uj is given by means of the cumulative chord-length method [12]. The optimization process is then divided into two steps. These two steps are carried out repeatedly until the average error eav is converged.
C.-Y. Lin et al./Computers
16
In the first step, a least-squares method is employed to solve for the control points L?,.The B-spline curve is expressed as a set of simultaneous equations, where unknown control points Bi as given in Eq. (1) are dependent variables and the measurement points I; are independent variables. Eq. (1) can be rewritten in a homogenous matrix form as
PaTPI = [~lTEfwBl~
(7)
where [Y ] is an m X 3 matrix containing the Cartesian coordinates of the measurement points, [N] is a matrix with m rows and p + 1 columns containing the basis functions, and [B] is a (p + 1) X 3 matrix containing the unknown coordinates of the control points. Given uj with respect to each measurement point, the Gauss-Jordan elimination [13] is employed to solve Eq. (7) where the maximum dimension used in the computer program is m X ( p + 1). In the second step, the parameter values uj are optimized by means of the objective function given by Eqs. (5) and (6) to minimize the average error e a”* provided that the control points Bi are given. Since each e, as indicated in Eq. (6) is affected by uj only, the optimization process can be formulated as one-dimensional nonlinear optimization problem where each uj is optimized individually by the Brent’s search method [131. The Brent’s search method repeats itself for all measurement points to evaluate all uj and it finally gives the minimum of eav.
4. Surface lofting algorithm 4.1. Generation
of the lofted sw$ace
Consider a set of parallel curves obtained from the above curve fitting algorithm, denoted {P,(u),P,(u), . . . ,P,
in
Industry34 (1997)73-85 where Bi / represents the ith control point with respect to’the /th curve. A B-spline surface Q(u,u) can be defined as: p+l
q+1
where C,,j are the control points with respect to the B-spline surface, (q + 1) is the number of control segments along the u direction and M, k are the basis functions in the u direction. If the parameter u is fixed, i.e., u = ua, then Q(u, u) can be reduced to an isoparametric curve, as given below:
i=
(10) 1
q+l
T(‘O) = C
ci,,Mj,k(
‘0).
(11)
j= I
By comparing Eqs. (8) and (lo), we find that Wi(uO) can be considered as the control points of the isoparametric curves Q(u,u,). Therefore, if the surface Q(u,u) is divided into /‘,, isoparametric curves and W,(u,) are designed in such a way that Wi(uO) = where u0 E [O 11 and /= 1,. . ,/,, then all Bi,ft isoparametric curves Q(u,u,) would match the bestfitted curves P/(U). The so called skinning or lofting is used for fitting a surface to known curves. The surface should ideally pass through these curves in a spatial topology sequence. It is noted that from the point of view of computer-aided design, a spline curve should generally be provided for guiding the generation of the lofted surface [5,6]. For the reverse engineering application, however, no spline curve is available, and, hence the available algorithms as indicated in the literature are invalid. Instead of using this approach, we propose another method by employing the curve fitting algorithm as described in Section 3 for fitting the control points C;,j to known points Bi,/ for each i. This can ensure that the created surface to be visually smooth. A detailed description of the algorithm for optimizing the control points Ci j is described below. Consider a set of curves P;(U), /= 1,. . ,/,, as given in Eq. (8). To clarify the explanation a case with four curves, each of which has four control points, is illustrated, as depicted in Fig. 1. Each set
C.-Y. Lin et al. /Computers in Industry 34 (19971 73-85
II
less than that of Bi,/, the curve Cj would become smoother. For this situation Ci does not pass through Bi,/, and, hence Eq. (12) must be revised as follows: Bi ./= Cd %,e> 7
(13)
where Ci(vi /> represents a particular point on the curve Ci wlkch is most close to the point B,,/. The overall surface accuracy can properly be controlled by optimizing the control points Ci,j, and, hence the curves Ci by means of the curve fitting algorithm as described in Section 3 4.2. Determination
of the parameter
v
The formulas for the lofted surface can be summarized as follows: p+l Q(
uTui,/)
=
q+1
C i=
C 1
[
j=
‘i,j”j,k( 1
Pfl
Fig. 1. A case with 4 curves, control points.
each of which is composed
of 4
=
C i=
‘i,/)
1 Ivl,k(
‘)
ci(vi,/)&,k(u)* 1 (14)
of the control points Bi,/ with fixed i is fitted by a curve, denoted C by means of the curve fitting algorithm as described in Section 3. A set of parametric values ui,/ corresponding to the points B,,/ on the curves Ci cm also be determined, i.e., B;.P= ‘i( 'i,f)
.
(12)
Consider the curve C, as shown in Fig. 1 as an example, where C, is fitted from the points B,,,, B 2.2, B2.3 and B2.4. The control points of C, are C 1 c22, C,, and Cz4, and the parametric values ,f;h respect ;o B,,, , b,,, , 4, and 5, are u2,1’ respectively. The number of the md ‘2,4’ ‘2.2’ ‘2,3 control points Ci,j does not have to be equal to that of the points Bi,/. If both are equal, i.e. j =/‘, the curve Ci would completely pass through B,,/. However this might result in the crease problem as the points Bi./ are the control points of the curves P,(U) which are fitted from the measurement data, rather than extracted from an ideal surface model. On the contrary, if the number of the control points Cj,j is
Note that the parametric values vi,/ in the above equation are discrete and are defined only on the points Ci(vi,/) which are most close to the points Bi,/. But, since the u and v parameters in a surface model should vary continuously, as evidenced in Eq. (lo), it is therefore necessary to define a continuous function for the v parameter so that any point on the surface Q(u,v) can be defined. Consider the case as shown in Fig. 2, where there are 4 X 4 discrete values vi,/ corresponding to 4 X 4 control points B,,/. The vi+/ values are listed in a rectangular array as follows:
‘i,f
e=1 e= 2 f= 3 e=4
i==l 0 0.3 0.6 1
i=2 0 0.32 0.58 1
i=3 0 0.28 0.6 1
i=4 0 0.3 0.62 1
This table indicates that when L’= 2 and 3, the vi./ values along the u direction are unequal. This means that the u parameter along any section curve in the u direction may vary. For example, the interpolation
Lin et
Computers
73-85
curve P2(u), is the ratio between Q,Qz and the total length, i.e., VI = 0.3 in this case. Similarly, the equivalent v value, us, for the third curve P,(u) is the ratio between Q, Q2 + Q2 Q3 and the total length, i.e., v: = 0.6. In summary, the equivalent values u; for the four curves P,, P2, P, and P4 as indicated in Fig. 3 are 0, 0.3, 0.6 and 1, respectively. Once v; are defined, the function fi(vi,/> is given as the following form: fi(vi,!)
=0
forall
i= l,...,
n+
1,if
v’=O, (16)
f.(ui,/)
= 1
forall
i=
l,...,
n+
1,if
uc=
1,
( 17) otherwise,
Fig. 2. The parameter values ui,/ are discretely distributed. An interpolation function must be added to estimate the values ui for the ikerpolated curve.
curve as depicted by 2 has four different the u direction. The u parameter face can be described
vi,;=fxVi,/)*
a dashed curve as shown in Fig. v values, v’, v*, v3 and v4 in + of any point on the lofted surby the following equation:
(15)
where the subscript i at vi,; represents that the v parameter with respect to a point on the surface varies not only in the u direction, but also in the u direction. The interpolation function h(vi,/) is used to provide a continuous variation of the v parameter over the entire lofted surface. Note that two constraints should be specified for the function L.(u~,~>. First, fi(ui,/) must be one-to-one mapping. Second, fi(vi,/> must at least be C’ continuous so that the smoothness of the lofted surface can be sustained. Fig. 3 illustrates the method for determining the function fi(~~,~). The centroid of each curve P,(u) is obtained first, denoted Q,. The equivalent v value at the centroid, denoted v>, for each curve is then determined based on the chord length method [12]. In Fig. 3, Q,, Q2, Q3 and Q4 represent the centroids of P,(u), PJu), P,(u) and P4(u), respectively. The -total length of the four centroids is Q, Q2 + Q2 Q3 + Q3Q4. The equivalent v value, vi, for the second
(
L;,,/-
u; “,,e+ 1 -
L$+, -
“i,L
v;
)I (18) ’
where vc represents the equivalent v value of the centroid of the interpolated curve and d is a con-
Fig. 3. The distribution
of the centroid for each curve.
C-Y. Lin et al./ Computers in Industry 34 (1997) 73-85
79
the Levenberg-Marquardt method [3]. Hereafter this approach is termed B-spline surface fitting. 5.1. Computer simulation
Fig. 4. The shape of the interpolation
function
fi.
stant, as indicated in Fig. 4. Eqs. (16)-(18) basically provide a linear interpolation for specifying the v value in the u direction. But, the intersection between two adjacent line segments is made smooth to ensure that the function fi(ui !) is C’ continuous. Consider the second curve C,’ as depicted in Fig. 2 as an example. The u value of any point on C, can be specified by the function f&~). The shape of f2(v2,e> is depicted in Fig. 4, where the horizontal axis represents the equivalent u value of the centroid of the curve P,(u) and the vertical axis represents the u value of the intersection point between C, and the interpolated curve. The parameter d determines the smoothness of the function fi(uz,e>. If d is too small, the smoothness of f2(u2,/) becomes worse. But, if d is too large, the deviation between ui,j and V , e when uc = v; may increase.
For the computer simulation, the data points used are generated by a simulation program in which a random number is added to each of the x, y and z components of the points with respect to an ideal CAD model to simulate a practical condition. All simulations are carried out on an IBM PC 486/33 MHz with 8MB RAM. In the first example, an isolated surface of a size 30 X 20 mm with control points 7 X 5 is used to demonstrate the feasibility of the proposed surfacelofting method. The number of the simulated measurement points is 40 X 40, with a maximum random value 0.01 mm. Fig. 5 depicts the distribution of the errors between the measurement points and the lofted surface, where the errors are amplified by 10 times. The control points used for both methods are 10 X 10. A quantitative comparison of these two approaches is indicated in Table 1. It indicates that for the surface lofting method, the average error before optimization is 3.2360 mm which has been reduced to 0.0204 mm and 0.0178 mm, respectively, in 10 and 15 iterations, meanwhile, for the B-spline surface fitting, the average error starts from 0.0395 mm and converges to 0.0166 mm and 0.0130 mm, respectively, in 10 and 15 iterations. Here the average error e av is defined as I1
m
n
(19) ejj = Iqj - P( uiuj)l,
(20)
5. Examples Several examples including computer simulation and experiments are presented in this section to illustrate the feasibility of the proposed surface-lofting approach. For the purpose of comparison, a surface optimization method [ 141 is employed also to fit the measurement data by a B-spline surface model as given in Eq. (911. This B-spline surface fitting method has been demonstrated to be computationally more efficient and less CPU memory required than
Fig. 5. Errors between the simulated measurement points and the lofted surface, where the scale factor for errors is 10.
80
C.-Y. Lin et al. / Computers in Industry 34 (1997) 73-85
Table 1 Results of optimization
for a sculptured
surface with 40 X 40 measurement
points; the control points used in the fitted surface are 10 X 10 B-spline surface fitting
Surface lofting
Iterations numbers
eav Before optimization 10 iterations 15 iterations
(mm)
Effective computing
3.2360 0.0204 0.0178
rate
0.67 1”
eau (mm)
Effective computing
0.0395 0.0166 0.0130
2.36 3.43
rate
* The CPU time required for tbis case is 1090 s and is acted as the standard. Effective computing rate = (CPU time required)/1090.
where m and IZ represent the number of measurement points along the n and y directions, respectively, eij is the shortest distance between the measurement point Yij and the lofted surface. Table 1 indicates that the average error for both methods can be reduced to the same level, with the B-spline surface fitting being slightly better than the surface lofting method. However, the effective computing rate, as defined in Table 1, for the surface lofting method (0.67 and 1 in 10 and 15 iterations, respectively) is significantly less than that of the B-spline surface fitting (2.36 and 3.43 in 10 and 15 iterations, respectively). This result indicates that the surface lofting method is computationally more efficient than the B-spline surface fitting. Table 2 lists the results of optimization for a various number of control points for the same case as shown in Table 1. This table indicates that the average error is reduced from 0.0779 to 0.0046 mm, respectively, when the control points are increased
from 6 X 6 to 11 X 11 for the surface lofting method, meanwhile, it is reduced from 0.0661 to 0.0041 mm for the B-spline surface fitting. However, it should be noted that the surface would oscillate when too many control points are used. Table 2 also indicates that the CPU time required for the B-spline surface fitting is more sensitive to the number of control points. For example, when the control points are increased from 6 X 6 to 11 X 11, the effective computing rate is only increased from 0.69 to 1.06 for the surface lofting method, but it is increased from 1.54 to 4.26 for the B-spline surface fitting. In the second example, an isolated surface of a size 30 X 20 mm with control points 7 X 5 is used. The number of the simulated measurement points is 60 X 60, with a maximum random value 0.02 mm. Table 3 lists the results of optimization for various number of patches. This table indicates that for the B-spline surface fitting the stack would be overflow when the measurement points exceed 60 X 30;
Table 2 Effects of control points on the average error eaV and the CPU time required for the surface as indicated Number of iterations: Control points
6X6 7x7 8X8 9x9 10 x 10 11 x 11 * Effective computing
in Table 1
15 Surface-lofting
method
e Gfore
eav after
(mm)
hd
3.8904 3.7248 3.5577 3.3942 3.2136 3.0786
0.0779 0.0445 0.0347 0.0131 0.0178 0.0046
B-spline surface fittingpoints Effective computing rate *
eav before optimization
0.69 0.78 0.87 0.98 1 1.06
0.1357 0.1415 0.0759 0.0390 0.0395 0.0083
rate is defined in Table 1.
(mm)
eav after optimization
0.0661 0.0402 0.0284 0.0121 0.0130 0.0041
(mm)
Effective computing
1.54 1.90 2.34 2.90 3.43 4.26
rate
*
C.-Y. Lin et al. / Computers in Industry 34 (1997) 73-85
whereas it is still possible for the surface lofting method to describe the above surface with one patch. Fig. 6(a) illustrates that 2 X 2 patches, each of which is a B-spline surface model, are used to describe the original shape. Two adjacent surfaces are then connected by means of a surface blending algorithm as given in Lai et al. [14]. Fig. 6(b) depicts the one-patch surface model by means of the proposed surface-lofting method. The average error for all three cases, i.e., 4, 2 and 1 patches, for the !surface lofting method is 0.0207 mm, 0.0193 mm a.nd 0.0203 mm, respectively, as indicated in Table 3; whereas for the B-spline surface fitting it is 0.0169 mm. The maximum difference among these values is only 0.0038 mm, which is very small and can be neglected. This same table also indicates that the CPU time required for the surface lofting method for all three cases is substantially less than the 15spline surface fitting. Although the effective computing rate for the surface lofting method is increased from 1.39 to 1.66 when the patch number decreases from 4 to 1, the case with less patch number is still suggested as the CAD data format would be simpler. The curve fitting algorithm used in the surface lofting method is ressentially one-dimensional, i.e., the fitting is carried out along one parametric direction. However, the B-spline surface fitting is two-di-
Table 3 Effects of the number of patches on the average error e, Number of iterations: Number of patches
81
Fig. 6. Fitted surface models for an example with 60x60 data points; (a) the B-spline surface fitting, where each patch has 5 X 4 control points; (b) the lofting surface, where the control points are 10X8.
and the CPU time required for a surface with 60 X 60 measurement
No.
Control points Effective computing
Number of data points
Surface-lofting method
B-spline surface fitting
eav (mm)
Effective computing
eav (mm)
rate *
rate *
4
1 2 3 4
5x4 5x4 5x4 5x4
30 x 30 30 x 30 30x30 30x30
0.32 0.36 0.35 0.36
0.0213 0.0265 0.0166 0.0168
0.76 0.78 0.85 0.85
0.0164 0.0210 0.0144 0.0152
2
Sum 1 2
10 x 4 10 x 4
60 x 30 60 x 30
1.39 0.85 0.81
0.0207 0.0221 0.0160
3.25 ** **
0.0169 ** *
Sum 1
10x 8
60X60
1.65 1.66
0.0193 0.0203
** **
1
* Effective computing * ’ Not executive.
points
15
rate is delined in Table 1.
l
l
l
* I
82
C.-Y. Lin et al./ Computers in Industry 34 (1997) 73-85
mensional, where the fitting is carried out along the u and u directions simultaneously. The parameter optimization processes for both algorithms are similar, i.e., Gauss-Jordan elimination and Brent’s search method. But, the maximum dimensions for the curve fitting is m X (p + I), as indicated in Eq. (7), whereas it is (mn) X [( p + 1Xq + l)] for the B-spline surface fitting, as indicated in Lai et al. [ 141. Accordingly, the B-spline surface fitting requires much more arithmetic operations and CPU memory, as evidenced from previous two examples. 5.2. Experiments A cartoon model with a size 80 X 90 mm is used to verify the feasibility of the proposed approach for complex shape modeling. A Keyence LC-2030 laser displacement sensor is mounted vertically on the z-axis of an x-y-z table. This probe has a resolution of 0.5 p,rn, a reference range of 51 mm and a measuring range of f8 mm. The digitizing process is undertaken by application of the raster scan mode. That is, the laser probe moves horizontally along one direction, say x, and then increases a step size along the y direction. It moves in a constant speed and takes the data periodically. A step size of 0.5 mm is used both along the x and y directions. Fig. 7 depicts the measurement data for one quarter of the cartoon model, where there are 40 lines and each has 185 points. This figure clearly indicates that the data points are scarce near the lower part of the plot, which is because a sudden fall appears in this region. Fig. 8 illustrates the curve fitting results for the above case, where each curve is fitted by 30
Fig. 8. Unexpected oscillation appears in the region with scarce data points, where each curve is fitted by 30 control points.
control points and two of the curves are drawn independently to illustrate the wave form. The region marked by a circle as indicated in Fig. 8 is supposed to be flat, but, unexpected oscillation appears due to the scarce distribution of the data points. Even the increase in the number of control points cannot completely eliminate such a phenomenon. A linear interpolation method is employed here to resolve the unexpected oscillation problem on the curves generated. The purpose of this method is to interpolate a number of points between any two adjacent measurement points so as to reduce the length of the line segment between two adjacent points. Consider that Q, and Q, are two adjacent measurement points, Pi is the ith interpolation point between Q, and Q, and dmin is a pre-defined tolerance. Then, Pi can be determined by means of the following equations: k=
IQ,-Q,l dmin
P,=e,+$(e,-e,),
Fig. 7. The measurement data for one quarter of the cartoon model, where there are 40 lines and each has 185 points.
where k tween Q, oscillation which is
k is an integer,
(21)
’
i=
l,...,k,
(22)
is the number of interpolation points beand Q,. Fig. 9(a) depicts the unexpected that occurs on the lower part of the curve generated by means of original measure-
C-Y. tin et ai. / Computers
in Indusrry 34 (1997) 73-8.5
83
ment points. However, after the linear interpolation is implemented, as jndicated in Fig. 9(b), the oscillation disappears. Fig. 10 illustrates that the unexpected oscillation near the lower part of the plot as indicated in Fig. 8 has become flat after the liner interpolation is implemented. Therefore, the proposed linear interpolation method can provide an easy manner for eliminating the oscillation problem due to too scarce data points on the part surface. Fig. 11 depicts the complete CAD model of the cartoon model by mltans of the surface lofting method when the measurement points are fitted along the u
I1 r
Control point
--The
A
Fig. 10. Unexpected oscillation as indicated in Fig. 8 disappears after the linear interpolation is implemented, where the control points are 15 X 25 and the number of iterations is 5.
\--Original
cwve fitted
measurement point
*-
_
The curve fitted
direction, where there are 185 lines and each has about 160 points. This CAD model is composed of three patches, each of which has 55 X 35 control points. The dark lines between two adjacent surface patches as indicated in Fig. 11 represent the cubic spline blending surface, used to connect adjacent surface patches together [ 141. The error between each measurement point and the fitted curves can be computed. In this case the average error between the measurement points and the fitted curves was observed to be less than 0.02 mm. It is noted that when the B-spline surface fitting is employed for this case, 5 X 5 patches must be used. Also, note that for the surface lofting method the number of data points on each scan line is allowed to be different, whereas it should be equal for the B-spline surface fitting. This makes the proposed method more applicable for practical applications.
6. Conclusion Interpolated point Original measurement point
(b) Fig. 9. (a) The fitted curve before the linear interpolation is implemented, (b) The fitted curve after the linear interpolation is implemented. The number of control points is 15 and the number of iterations is 5.
A surface lofting method is proposed in this work for the reverse engineering of complex shape. This approach includes a curve fitting algorithm and a surface lofting algorithm, where the former is used to fit for a series of section curves from the measurement data and the latter is used to fit for the lofted surface from the section curves. Tbe surface lofting
C.-Y. Lin et al./Computers
in Industry 34 (1997) 73-M
Fig. Il. The complete lofted surface model for the cartoon model, where 3 patches are used to describe the entire surface.
algorithm is essentially a curve fitting where each set of the control points, corresponding to the same index on each of the section curves, is fitted to a new curve. The control points of all the new curves form the control points of the desired surface. The above two fitting algorithms are essentially one-dimensional (curve fitting), whereas the B-spline surface fitting is two-dimensional. This makes the proposed method more effective in terms of computational efficiency and CPU memory required, as evidenced from the comparison study in previous section. In addition, the proposed method is allowed to have different number of measurement data on each scan line, which makes it more applicable for arbitrary shape modeling. Moreover, the linear interpolation algorithm proposed in this work can effectively eliminate the distortion of the lofted surface near the sharp edges and comers. The proposed approach can be extended for the reverse engineering of other surface models such as surface of revolution, sweep surface and ruled surface. This may enhance the capability of the proposed approach for composite surface modeling. Other related works such as boundary detection of the surface patches and surface integration are also worthy of studying.
Acknowledgements The authors would like to acknowledge the National Science Council, ROC, for its support through grant No. NSC-86-2212-E008-008. The authors would also like to acknowledge the reviewers for useful comments.
References [l] J. Hoschek, Spline approximation of offset curves, Computer Aided Geometric Design 5 (1988) 33-40. [2] D.F. Rogers, N.G. Fog, Constrained B-spline curve and surface fitting, Computer Aided Geometric Design 21 (10) (1989) 641-648. [3] B. Sarkar, C.H. Menq, Smooth-surface approximation and reverse engineering, Computer Aided Geometric Design 23 (9) (1991) 623-628. [4] B. Sarkar, C.H. Menq, Parameter optimization in approximating curves and surfaces to measurement data, Computer Aided Geometric Design 8 (1991) 267-290. [5] W. Tiller, Rational B-splines for curve and surface representation, IEEE Comput. Graphics and Applications 3 (10) (1983) 61-69. [6] C.D. Woodward, Skinning techniques for iterative B-spline interpolation, Computer Aided Geometric Design 20 (8) (1988) 441-451.
C.-Y. Lin et al./ Computers in Industry 34 (1997) 73-85 surfaces as [71F. Lin, W.T. Hewitt, Expressing Goons-Gordon NURBS, Computer Aided Geometric Design 26 (2) (1994) 145-155. [81D.J. Filip. T.W. Ball, Procedurally representing lofted surfaces, IEEE Comput. Graphics and Applications (1989) 2733. [91E. Cohen, T. Lynche, R.F. Riesenfeld, Discrete B-spline and sub-division techniques in computer aided geometric design and computer graphics, Comput. Graphics and Image Processing 14 (1980) 137-l 11. I101W. Bohem, H. Prautzsch, The insertion algorithm, Computer Aided Geometric Design 17 (1) (1985) 58-59. Algorithm for [illE. Cohen, T. Lyche, L.L. Schumarker, degree-raising of splines, ACM Trans. Graph. 4 (3) (1985) 171-181. [I21D.F. Rogers, J.A. .4dams, Mathematical Elements for Computer Graphics, McGraw-Hill Publishing Company, New York, 1990. D31 G.V. Reklaitis, A. Ravindran, K.M. Ragsdell, Engineering Optimization - Methods and Applications, John Wiley and Sons, New York, 1984. [141J.-Y. Lai, C.-Y. Lin, Reverse engineering of composite sculp tured surfaces, Int. Journal of Advanced Manufacturing and Technology (1996) forthcoming. Cbih-Young Lin is a Ph.D. student in the Department of Mechanical Engineering at the National Central University, Taiwan, ROC. His research has focused on reverse engineering, quality control and management information system.
85 Chung-Shan Liou is a master student in the Department of Mechanical Engineering at the National Central University, Taiwan, ROC. His research has focused on CAD/CAM and reverse engineering.
Jiing-Yih Lai is a Professor in the Department of Mechanical Engineering at the National Central University, Taiwan, ROC. He obtained his M.Sc. and Ph.D. in Mechanical Engineering from the Ohio State University, Columbus, OH, USA. His research work has centered on CAD/CAM, reverse engineering and coordinate measuring metrology.
A surface-lofting
approach for smooth-surface reconstruction from 3D measurement data
Chih-Young Lin, Chung-Shan Liou, Jiing-Yih Lai Mechanical Engineering Department,
*
National Central Uniuersity, Chung-Li 320, Taiwan, ROC
Received 18 September 1995; revised 1 May 1996; accepted 27 August 1996
Abstract A reverse engineering method for sculptured surface modeling is proposed in this work. It includes a curve fitting algorithm and a surface lofting algorithm, where the former is used to fit each row of measurement data to obtain a B-spline curve and the latter ifs used to fit the section curves into a lofted surface. The above two fitting algorithms are essentially one-dimensional, whereas the commonly used B-spline surface fitting is two-dimensional. This makes the proposed method more effective in temis of computational efficiency and CPU memory required. Both computer simulation and experimental results indicate that I:he proposed method can efficiently be used for the reverse engineering of complex object shape. 0 1997 Elsevier Science B.V. Keywords: Reverse enguieering: Surface lofting; Curve fitting
1. Introduction Reverse engineering is a technology to establish CAD models from samples, prototypes, molds or manufactured parts by digitization. It can be used as a design tool for extracting the design concept of an existing product by analyzing its CAD model. The integration of rev’erse engineering and the CAM technology can even replace conventional copy milling method, commonly used to fabricate products from clays or wood models. There have been considerable studies with respect to fitting sculptured surfaces to known 3D measurement data [l-4]. Many problems still exist due to the complexity and variety of the shape of practical objects. Some of these are
* Corresponding author.
such as computational problem, the difficulty in processing enormous amount of data, the accuracy problem in created surface model, the difficulty for complex shape modeling, etc. The objective of this work is to present a surface lofting approach for the reverse engineering of complex shape. The problem of fitting parametric surfaces to known measurement data has been studied extensively. Numerous investigators have focused on the development of various surface optimization schemes for accurate modeling of sculptured surfaces. Some of the numerical methods proposed are such as the Newton-like approach [l], the first-order Taylor series expansion method [2] and the LevenbergMarquardt method [3]. These methods can generally yield results with comparable surface accuracy for an isolated sculptured surface. However, the Newton-
0166-3615/97/$17.00 0 1997 Elsevier Science B.V. All rights reserved. PII SOl66-3615(96)00082-6
74
C.-Y. Lin et al. / Computers in Industry 34 (19971 73-85
like approach [l] and the first-order Taylor series expansion method [2] are computationally less efficient than the Levenberg-Marquardt method 131. But, the Levenberg-Marquardt method requires a huge CPU memory, which limits its application to small amount of measurement data. The main drawback of the above surface fitting methods are that when the shape of the object becomes complex, none of them can yield effective CAD models. The maximum dimension used in the Levenberg-Marquardt method is (mn) X (mm), where m and n are the number of measurement points along the x and y coordinate axes, respectively. Both the computer time and CPU memory would increase substantially when the number of measurement points increases. Sarkar et al. [4] extended the Levenberg-Marquardt method by introducing a boundary detection scheme to divide the compound surface into several regions, each of which was fitted by a surface patch individually. However, how adjacent patches were connected was not addressed. The robustness issue of the proposed boundary detection scheme was not investigated also. A computer-aided design technique called ‘skinning’ or ‘lofting’ could be employed for the reverse engineering applications. This technique interpolates a family of curves defined in one of the surface’s parametric directions to form relatively complex surfaces. Tiller [5] and Woodward [6] gave a detailed description of constructing a lofted surface with Bsplines as well as the necessary background about interpolation. Lin et al. [7] presented a B-spline Boolean-sum approach for surface interpolation. A bivariate B-spline surface was constructed from an adequate superposition of univariate B-spline interpolation by the Boolean sum operation. The procedure was the same as that of Coons-Gordon surfaces. A possible way of extending the above approach to NURBS is also discussed for the approximation of Coons-Gordon-type surfaces. Filip et al. [8] pointed out several problems with using tensor-product B-spline surfaces in lofting, such as difficulty in lofting non-planar curves, how to determine appropriate knot vectors and the ‘crease’ problem. The crease problem arises when the surface normal does not change continuously, which results in a surface that is not visually smooth. They proposed a knot vector function, instead of a tensor-
product form where a single u knot vector is specified for every u value, such that the problem of lofting non-planar curves can be resolved. A procedure was also described for reparametering a G’ curve into a C’ curve so that the crease problem can be eliminated. An important characteristic of a set of curves capable to be lofted is that they should be comparable in degrees, number of spans and knot vector form. Cohen [9] proposed an Oslo algorithm for adding new knots. Associated algorithm for attaining new control points was provided also. Boehm [lo] developed another algorithm for adding knots and calculating new control points. Such an algorithm is actually a special case of the Oslo algorithm, but is computationally more efficient. Cohen et al. [I I] proposed a degree-raising approach for adding flexibility to a set of B-spline curves.
2. Problem statement The primary limitation in the available surface optimization methods [l-4] is that they cannot resolve the problem of complex shape modeling. When the measurement data increases, the substantial increase in computational time and CPU memory makes these methods ineffective. Moreover, a single isolated surface is often not good enough to describe the shape of the object surface. It is generally necessary to divide the object surface into a number of small patches such that the surface optimization algorithm can be applied. The objective of this work is to develop a method incorporating a curve fitting algorithm and a surface lofting algorithm for the reverse engineering of complex physical parts. The curve fitting algorithm is employed first to fit each row of measurement data, which yields a set of parallel B-spline curves. A least-squares algorithm is employed to optimize the control points and parameter values for the B-spline curves, so as to minimize the average deviation between each row of measurement data and the curve fitted. The surface lofting algorithm is then employed for interpolating the above parallel curves to construct a lofted surface. The proposed surface-lofting algorithm is different from those described in the literature [5-81. The
C.-Y. Lin et al. /Computers in Industry 34 (1997) 73-85
main differences are explained below. Firstly, for the CAD application, fhe describing curves are designed in the CAD system and the lofted surface is expected to pass all describing curves. For the reverse engineering application, however, the created surface would distort if all section curves are to be passed. This is because the available surface-lofting algorithms are not suitable for curves which are densely distributed. A fittin,g algorithm is developed to fit the section curves into a smooth surface. Secondly, the surface accuracy is an important issue in surface reconstruction. The lofted surface must be able to describe all points on the part surface. A reparameterization process :is employed to match this goal. Finally, sharp edges and comers appear frequently in the reverse engineering problem, which may lead to unexpected oscillation on the created surface. A detail discussion o-F the above issues will be given and methods will be provided to solve these problems.
of U. The basis functions Ni,k(~) are related to the knot vector U, and are defined as [12]: 1 N,,I =
ifx;
i 0
N, (Uj =
fitting algorithm
The curve fitting problem has been extensively studied and various computational algorithms are available in the literature [1,2,4]. However, for the completeness of thi;s paper, the algorithm used in this work is briefly described. Each array of measurement points (along the x direction) is modeled by a B-spiine curve. The measurement points are densely distributed such that the cross sectional profile of the part surface can truly be described. An opeln cubic B-spline curve is used, which can handle b’oth smoothness and local control property for each array of measurement data to be fitted. The position vector P(U) of a B-spline curve along the parameter u can be expressed as P+l
2
1,
(1)
where B, are the position vectors of the control points, p + 1 is the number of the control segments, P&(U) are normalized B-spline basis functions and urnin and s,,, are ,the lower and the upper margins
(2)
otherwise, (“-xi>N,,k-l(u)
r,k
x;+!c-1
+
(Xifk
--xi
-
wi+wG4
xi+
k -xi+
3
(3)
1
where xi are the knot values of the knot vector X. The normalized form of a kth order open uniform knot vector X is given as [12]: X =
[ x1 x2
. . . xkXk+l
xi = 0, 1 I ilk;
... x,+1 x,+2 ... x,+,‘+l
xi =
k+l
3. Curve
15
1f (4)
i-k p+2-k’
xi=1,p+2~i
The curve fitting algorithm is employed to optimize the position of the control points Bi and the parameter value u, as given in Eq. (11, corresponding to each measurement point. The other parameters, e.g., knot vectors and the number of control segments, in the B-spline model are determined beforehand. The average error eav is defined as the average of the sum of the squares of the fitting errors at each data point, i.e. e a” ZZ
;)!,
Cej)’
’
(5)
Je,=Iq-P(uj)I,
(6)
where m is the number of the measurement points, ej represents the shortest distance between the jth measurement point 5 and the curve fitted. The parameter optimization problem is formulated as the minimization of the error function eaV which yields an optimum set of solutions for Bi and uj. An initial estimation of uj is given by means of the cumulative chord-length method [12]. The optimization process is then divided into two steps. These two steps are carried out repeatedly until the average error eav is converged.
C.-Y. Lin et al./Computers
16
In the first step, a least-squares method is employed to solve for the control points L?,.The B-spline curve is expressed as a set of simultaneous equations, where unknown control points Bi as given in Eq. (1) are dependent variables and the measurement points I; are independent variables. Eq. (1) can be rewritten in a homogenous matrix form as
PaTPI = [~lTEfwBl~
(7)
where [Y ] is an m X 3 matrix containing the Cartesian coordinates of the measurement points, [N] is a matrix with m rows and p + 1 columns containing the basis functions, and [B] is a (p + 1) X 3 matrix containing the unknown coordinates of the control points. Given uj with respect to each measurement point, the Gauss-Jordan elimination [13] is employed to solve Eq. (7) where the maximum dimension used in the computer program is m X ( p + 1). In the second step, the parameter values uj are optimized by means of the objective function given by Eqs. (5) and (6) to minimize the average error e a”* provided that the control points Bi are given. Since each e, as indicated in Eq. (6) is affected by uj only, the optimization process can be formulated as one-dimensional nonlinear optimization problem where each uj is optimized individually by the Brent’s search method [131. The Brent’s search method repeats itself for all measurement points to evaluate all uj and it finally gives the minimum of eav.
4. Surface lofting algorithm 4.1. Generation
of the lofted sw$ace
Consider a set of parallel curves obtained from the above curve fitting algorithm, denoted {P,(u),P,(u), . . . ,P,
in
Industry34 (1997)73-85 where Bi / represents the ith control point with respect to’the /th curve. A B-spline surface Q(u,u) can be defined as: p+l
q+1
where C,,j are the control points with respect to the B-spline surface, (q + 1) is the number of control segments along the u direction and M, k are the basis functions in the u direction. If the parameter u is fixed, i.e., u = ua, then Q(u, u) can be reduced to an isoparametric curve, as given below:
i=
(10) 1
q+l
T(‘O) = C
ci,,Mj,k(
‘0).
(11)
j= I
By comparing Eqs. (8) and (lo), we find that Wi(uO) can be considered as the control points of the isoparametric curves Q(u,u,). Therefore, if the surface Q(u,u) is divided into /‘,, isoparametric curves and W,(u,) are designed in such a way that Wi(uO) = where u0 E [O 11 and /= 1,. . ,/,, then all Bi,ft isoparametric curves Q(u,u,) would match the bestfitted curves P/(U). The so called skinning or lofting is used for fitting a surface to known curves. The surface should ideally pass through these curves in a spatial topology sequence. It is noted that from the point of view of computer-aided design, a spline curve should generally be provided for guiding the generation of the lofted surface [5,6]. For the reverse engineering application, however, no spline curve is available, and, hence the available algorithms as indicated in the literature are invalid. Instead of using this approach, we propose another method by employing the curve fitting algorithm as described in Section 3 for fitting the control points C;,j to known points Bi,/ for each i. This can ensure that the created surface to be visually smooth. A detailed description of the algorithm for optimizing the control points Ci j is described below. Consider a set of curves P;(U), /= 1,. . ,/,, as given in Eq. (8). To clarify the explanation a case with four curves, each of which has four control points, is illustrated, as depicted in Fig. 1. Each set
C.-Y. Lin et al. /Computers in Industry 34 (19971 73-85
II
less than that of Bi,/, the curve Cj would become smoother. For this situation Ci does not pass through Bi,/, and, hence Eq. (12) must be revised as follows: Bi ./= Cd %,e> 7
(13)
where Ci(vi /> represents a particular point on the curve Ci wlkch is most close to the point B,,/. The overall surface accuracy can properly be controlled by optimizing the control points Ci,j, and, hence the curves Ci by means of the curve fitting algorithm as described in Section 3 4.2. Determination
of the parameter
v
The formulas for the lofted surface can be summarized as follows: p+l Q(
uTui,/)
=
q+1
C i=
C 1
[
j=
‘i,j”j,k( 1
Pfl
Fig. 1. A case with 4 curves, control points.
each of which is composed
of 4
=
C i=
‘i,/)
1 Ivl,k(
‘)
ci(vi,/)&,k(u)* 1 (14)
of the control points Bi,/ with fixed i is fitted by a curve, denoted C by means of the curve fitting algorithm as described in Section 3. A set of parametric values ui,/ corresponding to the points B,,/ on the curves Ci cm also be determined, i.e., B;.P= ‘i( 'i,f)
.
(12)
Consider the curve C, as shown in Fig. 1 as an example, where C, is fitted from the points B,,,, B 2.2, B2.3 and B2.4. The control points of C, are C 1 c22, C,, and Cz4, and the parametric values ,f;h respect ;o B,,, , b,,, , 4, and 5, are u2,1’ respectively. The number of the md ‘2,4’ ‘2.2’ ‘2,3 control points Ci,j does not have to be equal to that of the points Bi,/. If both are equal, i.e. j =/‘, the curve Ci would completely pass through B,,/. However this might result in the crease problem as the points Bi./ are the control points of the curves P,(U) which are fitted from the measurement data, rather than extracted from an ideal surface model. On the contrary, if the number of the control points Cj,j is
Note that the parametric values vi,/ in the above equation are discrete and are defined only on the points Ci(vi,/) which are most close to the points Bi,/. But, since the u and v parameters in a surface model should vary continuously, as evidenced in Eq. (lo), it is therefore necessary to define a continuous function for the v parameter so that any point on the surface Q(u,v) can be defined. Consider the case as shown in Fig. 2, where there are 4 X 4 discrete values vi,/ corresponding to 4 X 4 control points B,,/. The vi+/ values are listed in a rectangular array as follows:
‘i,f
e=1 e= 2 f= 3 e=4
i==l 0 0.3 0.6 1
i=2 0 0.32 0.58 1
i=3 0 0.28 0.6 1
i=4 0 0.3 0.62 1
This table indicates that when L’= 2 and 3, the vi./ values along the u direction are unequal. This means that the u parameter along any section curve in the u direction may vary. For example, the interpolation
Lin et
Computers
73-85
curve P2(u), is the ratio between Q,Qz and the total length, i.e., VI = 0.3 in this case. Similarly, the equivalent v value, us, for the third curve P,(u) is the ratio between Q, Q2 + Q2 Q3 and the total length, i.e., v: = 0.6. In summary, the equivalent values u; for the four curves P,, P2, P, and P4 as indicated in Fig. 3 are 0, 0.3, 0.6 and 1, respectively. Once v; are defined, the function fi(vi,/> is given as the following form: fi(vi,!)
=0
forall
i= l,...,
n+
1,if
v’=O, (16)
f.(ui,/)
= 1
forall
i=
l,...,
n+
1,if
uc=
1,
( 17) otherwise,
Fig. 2. The parameter values ui,/ are discretely distributed. An interpolation function must be added to estimate the values ui for the ikerpolated curve.
curve as depicted by 2 has four different the u direction. The u parameter face can be described
vi,;=fxVi,/)*
a dashed curve as shown in Fig. v values, v’, v*, v3 and v4 in + of any point on the lofted surby the following equation:
(15)
where the subscript i at vi,; represents that the v parameter with respect to a point on the surface varies not only in the u direction, but also in the u direction. The interpolation function h(vi,/) is used to provide a continuous variation of the v parameter over the entire lofted surface. Note that two constraints should be specified for the function L.(u~,~>. First, fi(ui,/) must be one-to-one mapping. Second, fi(vi,/> must at least be C’ continuous so that the smoothness of the lofted surface can be sustained. Fig. 3 illustrates the method for determining the function fi(~~,~). The centroid of each curve P,(u) is obtained first, denoted Q,. The equivalent v value at the centroid, denoted v>, for each curve is then determined based on the chord length method [12]. In Fig. 3, Q,, Q2, Q3 and Q4 represent the centroids of P,(u), PJu), P,(u) and P4(u), respectively. The -total length of the four centroids is Q, Q2 + Q2 Q3 + Q3Q4. The equivalent v value, vi, for the second
(
L;,,/-
u; “,,e+ 1 -
L$+, -
“i,L
v;
)I (18) ’
where vc represents the equivalent v value of the centroid of the interpolated curve and d is a con-
Fig. 3. The distribution
of the centroid for each curve.
C-Y. Lin et al./ Computers in Industry 34 (1997) 73-85
79
the Levenberg-Marquardt method [3]. Hereafter this approach is termed B-spline surface fitting. 5.1. Computer simulation
Fig. 4. The shape of the interpolation
function
fi.
stant, as indicated in Fig. 4. Eqs. (16)-(18) basically provide a linear interpolation for specifying the v value in the u direction. But, the intersection between two adjacent line segments is made smooth to ensure that the function fi(ui !) is C’ continuous. Consider the second curve C,’ as depicted in Fig. 2 as an example. The u value of any point on C, can be specified by the function f&~). The shape of f2(v2,e> is depicted in Fig. 4, where the horizontal axis represents the equivalent u value of the centroid of the curve P,(u) and the vertical axis represents the u value of the intersection point between C, and the interpolated curve. The parameter d determines the smoothness of the function fi(uz,e>. If d is too small, the smoothness of f2(u2,/) becomes worse. But, if d is too large, the deviation between ui,j and V , e when uc = v; may increase.
For the computer simulation, the data points used are generated by a simulation program in which a random number is added to each of the x, y and z components of the points with respect to an ideal CAD model to simulate a practical condition. All simulations are carried out on an IBM PC 486/33 MHz with 8MB RAM. In the first example, an isolated surface of a size 30 X 20 mm with control points 7 X 5 is used to demonstrate the feasibility of the proposed surfacelofting method. The number of the simulated measurement points is 40 X 40, with a maximum random value 0.01 mm. Fig. 5 depicts the distribution of the errors between the measurement points and the lofted surface, where the errors are amplified by 10 times. The control points used for both methods are 10 X 10. A quantitative comparison of these two approaches is indicated in Table 1. It indicates that for the surface lofting method, the average error before optimization is 3.2360 mm which has been reduced to 0.0204 mm and 0.0178 mm, respectively, in 10 and 15 iterations, meanwhile, for the B-spline surface fitting, the average error starts from 0.0395 mm and converges to 0.0166 mm and 0.0130 mm, respectively, in 10 and 15 iterations. Here the average error e av is defined as I1
m
n
(19) ejj = Iqj - P( uiuj)l,
(20)
5. Examples Several examples including computer simulation and experiments are presented in this section to illustrate the feasibility of the proposed surface-lofting approach. For the purpose of comparison, a surface optimization method [ 141 is employed also to fit the measurement data by a B-spline surface model as given in Eq. (911. This B-spline surface fitting method has been demonstrated to be computationally more efficient and less CPU memory required than
Fig. 5. Errors between the simulated measurement points and the lofted surface, where the scale factor for errors is 10.
80
C.-Y. Lin et al. / Computers in Industry 34 (1997) 73-85
Table 1 Results of optimization
for a sculptured
surface with 40 X 40 measurement
points; the control points used in the fitted surface are 10 X 10 B-spline surface fitting
Surface lofting
Iterations numbers
eav Before optimization 10 iterations 15 iterations
(mm)
Effective computing
3.2360 0.0204 0.0178
rate
0.67 1”
eau (mm)
Effective computing
0.0395 0.0166 0.0130
2.36 3.43
rate
* The CPU time required for tbis case is 1090 s and is acted as the standard. Effective computing rate = (CPU time required)/1090.
where m and IZ represent the number of measurement points along the n and y directions, respectively, eij is the shortest distance between the measurement point Yij and the lofted surface. Table 1 indicates that the average error for both methods can be reduced to the same level, with the B-spline surface fitting being slightly better than the surface lofting method. However, the effective computing rate, as defined in Table 1, for the surface lofting method (0.67 and 1 in 10 and 15 iterations, respectively) is significantly less than that of the B-spline surface fitting (2.36 and 3.43 in 10 and 15 iterations, respectively). This result indicates that the surface lofting method is computationally more efficient than the B-spline surface fitting. Table 2 lists the results of optimization for a various number of control points for the same case as shown in Table 1. This table indicates that the average error is reduced from 0.0779 to 0.0046 mm, respectively, when the control points are increased
from 6 X 6 to 11 X 11 for the surface lofting method, meanwhile, it is reduced from 0.0661 to 0.0041 mm for the B-spline surface fitting. However, it should be noted that the surface would oscillate when too many control points are used. Table 2 also indicates that the CPU time required for the B-spline surface fitting is more sensitive to the number of control points. For example, when the control points are increased from 6 X 6 to 11 X 11, the effective computing rate is only increased from 0.69 to 1.06 for the surface lofting method, but it is increased from 1.54 to 4.26 for the B-spline surface fitting. In the second example, an isolated surface of a size 30 X 20 mm with control points 7 X 5 is used. The number of the simulated measurement points is 60 X 60, with a maximum random value 0.02 mm. Table 3 lists the results of optimization for various number of patches. This table indicates that for the B-spline surface fitting the stack would be overflow when the measurement points exceed 60 X 30;
Table 2 Effects of control points on the average error eaV and the CPU time required for the surface as indicated Number of iterations: Control points
6X6 7x7 8X8 9x9 10 x 10 11 x 11 * Effective computing
in Table 1
15 Surface-lofting
method
e Gfore
eav after
(mm)
hd
3.8904 3.7248 3.5577 3.3942 3.2136 3.0786
0.0779 0.0445 0.0347 0.0131 0.0178 0.0046
B-spline surface fittingpoints Effective computing rate *
eav before optimization
0.69 0.78 0.87 0.98 1 1.06
0.1357 0.1415 0.0759 0.0390 0.0395 0.0083
rate is defined in Table 1.
(mm)
eav after optimization
0.0661 0.0402 0.0284 0.0121 0.0130 0.0041
(mm)
Effective computing
1.54 1.90 2.34 2.90 3.43 4.26
rate
*
C.-Y. Lin et al. / Computers in Industry 34 (1997) 73-85
whereas it is still possible for the surface lofting method to describe the above surface with one patch. Fig. 6(a) illustrates that 2 X 2 patches, each of which is a B-spline surface model, are used to describe the original shape. Two adjacent surfaces are then connected by means of a surface blending algorithm as given in Lai et al. [14]. Fig. 6(b) depicts the one-patch surface model by means of the proposed surface-lofting method. The average error for all three cases, i.e., 4, 2 and 1 patches, for the !surface lofting method is 0.0207 mm, 0.0193 mm a.nd 0.0203 mm, respectively, as indicated in Table 3; whereas for the B-spline surface fitting it is 0.0169 mm. The maximum difference among these values is only 0.0038 mm, which is very small and can be neglected. This same table also indicates that the CPU time required for the surface lofting method for all three cases is substantially less than the 15spline surface fitting. Although the effective computing rate for the surface lofting method is increased from 1.39 to 1.66 when the patch number decreases from 4 to 1, the case with less patch number is still suggested as the CAD data format would be simpler. The curve fitting algorithm used in the surface lofting method is ressentially one-dimensional, i.e., the fitting is carried out along one parametric direction. However, the B-spline surface fitting is two-di-
Table 3 Effects of the number of patches on the average error e, Number of iterations: Number of patches
81
Fig. 6. Fitted surface models for an example with 60x60 data points; (a) the B-spline surface fitting, where each patch has 5 X 4 control points; (b) the lofting surface, where the control points are 10X8.
and the CPU time required for a surface with 60 X 60 measurement
No.
Control points Effective computing
Number of data points
Surface-lofting method
B-spline surface fitting
eav (mm)
Effective computing
eav (mm)
rate *
rate *
4
1 2 3 4
5x4 5x4 5x4 5x4
30 x 30 30 x 30 30x30 30x30
0.32 0.36 0.35 0.36
0.0213 0.0265 0.0166 0.0168
0.76 0.78 0.85 0.85
0.0164 0.0210 0.0144 0.0152
2
Sum 1 2
10 x 4 10 x 4
60 x 30 60 x 30
1.39 0.85 0.81
0.0207 0.0221 0.0160
3.25 ** **
0.0169 ** *
Sum 1
10x 8
60X60
1.65 1.66
0.0193 0.0203
** **
1
* Effective computing * ’ Not executive.
points
15
rate is delined in Table 1.
l
l
l
* I
82
C.-Y. Lin et al./ Computers in Industry 34 (1997) 73-85
mensional, where the fitting is carried out along the u and u directions simultaneously. The parameter optimization processes for both algorithms are similar, i.e., Gauss-Jordan elimination and Brent’s search method. But, the maximum dimensions for the curve fitting is m X (p + I), as indicated in Eq. (7), whereas it is (mn) X [( p + 1Xq + l)] for the B-spline surface fitting, as indicated in Lai et al. [ 141. Accordingly, the B-spline surface fitting requires much more arithmetic operations and CPU memory, as evidenced from previous two examples. 5.2. Experiments A cartoon model with a size 80 X 90 mm is used to verify the feasibility of the proposed approach for complex shape modeling. A Keyence LC-2030 laser displacement sensor is mounted vertically on the z-axis of an x-y-z table. This probe has a resolution of 0.5 p,rn, a reference range of 51 mm and a measuring range of f8 mm. The digitizing process is undertaken by application of the raster scan mode. That is, the laser probe moves horizontally along one direction, say x, and then increases a step size along the y direction. It moves in a constant speed and takes the data periodically. A step size of 0.5 mm is used both along the x and y directions. Fig. 7 depicts the measurement data for one quarter of the cartoon model, where there are 40 lines and each has 185 points. This figure clearly indicates that the data points are scarce near the lower part of the plot, which is because a sudden fall appears in this region. Fig. 8 illustrates the curve fitting results for the above case, where each curve is fitted by 30
Fig. 8. Unexpected oscillation appears in the region with scarce data points, where each curve is fitted by 30 control points.
control points and two of the curves are drawn independently to illustrate the wave form. The region marked by a circle as indicated in Fig. 8 is supposed to be flat, but, unexpected oscillation appears due to the scarce distribution of the data points. Even the increase in the number of control points cannot completely eliminate such a phenomenon. A linear interpolation method is employed here to resolve the unexpected oscillation problem on the curves generated. The purpose of this method is to interpolate a number of points between any two adjacent measurement points so as to reduce the length of the line segment between two adjacent points. Consider that Q, and Q, are two adjacent measurement points, Pi is the ith interpolation point between Q, and Q, and dmin is a pre-defined tolerance. Then, Pi can be determined by means of the following equations: k=
IQ,-Q,l dmin
P,=e,+$(e,-e,),
Fig. 7. The measurement data for one quarter of the cartoon model, where there are 40 lines and each has 185 points.
where k tween Q, oscillation which is
k is an integer,
(21)
’
i=
l,...,k,
(22)
is the number of interpolation points beand Q,. Fig. 9(a) depicts the unexpected that occurs on the lower part of the curve generated by means of original measure-
C-Y. tin et ai. / Computers
in Indusrry 34 (1997) 73-8.5
83
ment points. However, after the linear interpolation is implemented, as jndicated in Fig. 9(b), the oscillation disappears. Fig. 10 illustrates that the unexpected oscillation near the lower part of the plot as indicated in Fig. 8 has become flat after the liner interpolation is implemented. Therefore, the proposed linear interpolation method can provide an easy manner for eliminating the oscillation problem due to too scarce data points on the part surface. Fig. 11 depicts the complete CAD model of the cartoon model by mltans of the surface lofting method when the measurement points are fitted along the u
I1 r
Control point
--The
A
Fig. 10. Unexpected oscillation as indicated in Fig. 8 disappears after the linear interpolation is implemented, where the control points are 15 X 25 and the number of iterations is 5.
\--Original
cwve fitted
measurement point
*-
_
The curve fitted
direction, where there are 185 lines and each has about 160 points. This CAD model is composed of three patches, each of which has 55 X 35 control points. The dark lines between two adjacent surface patches as indicated in Fig. 11 represent the cubic spline blending surface, used to connect adjacent surface patches together [ 141. The error between each measurement point and the fitted curves can be computed. In this case the average error between the measurement points and the fitted curves was observed to be less than 0.02 mm. It is noted that when the B-spline surface fitting is employed for this case, 5 X 5 patches must be used. Also, note that for the surface lofting method the number of data points on each scan line is allowed to be different, whereas it should be equal for the B-spline surface fitting. This makes the proposed method more applicable for practical applications.
6. Conclusion Interpolated point Original measurement point
(b) Fig. 9. (a) The fitted curve before the linear interpolation is implemented, (b) The fitted curve after the linear interpolation is implemented. The number of control points is 15 and the number of iterations is 5.
A surface lofting method is proposed in this work for the reverse engineering of complex shape. This approach includes a curve fitting algorithm and a surface lofting algorithm, where the former is used to fit for a series of section curves from the measurement data and the latter is used to fit for the lofted surface from the section curves. Tbe surface lofting
C.-Y. Lin et al./Computers
in Industry 34 (1997) 73-M
Fig. Il. The complete lofted surface model for the cartoon model, where 3 patches are used to describe the entire surface.
algorithm is essentially a curve fitting where each set of the control points, corresponding to the same index on each of the section curves, is fitted to a new curve. The control points of all the new curves form the control points of the desired surface. The above two fitting algorithms are essentially one-dimensional (curve fitting), whereas the B-spline surface fitting is two-dimensional. This makes the proposed method more effective in terms of computational efficiency and CPU memory required, as evidenced from the comparison study in previous section. In addition, the proposed method is allowed to have different number of measurement data on each scan line, which makes it more applicable for arbitrary shape modeling. Moreover, the linear interpolation algorithm proposed in this work can effectively eliminate the distortion of the lofted surface near the sharp edges and comers. The proposed approach can be extended for the reverse engineering of other surface models such as surface of revolution, sweep surface and ruled surface. This may enhance the capability of the proposed approach for composite surface modeling. Other related works such as boundary detection of the surface patches and surface integration are also worthy of studying.
Acknowledgements The authors would like to acknowledge the National Science Council, ROC, for its support through grant No. NSC-86-2212-E008-008. The authors would also like to acknowledge the reviewers for useful comments.
References [l] J. Hoschek, Spline approximation of offset curves, Computer Aided Geometric Design 5 (1988) 33-40. [2] D.F. Rogers, N.G. Fog, Constrained B-spline curve and surface fitting, Computer Aided Geometric Design 21 (10) (1989) 641-648. [3] B. Sarkar, C.H. Menq, Smooth-surface approximation and reverse engineering, Computer Aided Geometric Design 23 (9) (1991) 623-628. [4] B. Sarkar, C.H. Menq, Parameter optimization in approximating curves and surfaces to measurement data, Computer Aided Geometric Design 8 (1991) 267-290. [5] W. Tiller, Rational B-splines for curve and surface representation, IEEE Comput. Graphics and Applications 3 (10) (1983) 61-69. [6] C.D. Woodward, Skinning techniques for iterative B-spline interpolation, Computer Aided Geometric Design 20 (8) (1988) 441-451.
C.-Y. Lin et al./ Computers in Industry 34 (1997) 73-85 surfaces as [71F. Lin, W.T. Hewitt, Expressing Goons-Gordon NURBS, Computer Aided Geometric Design 26 (2) (1994) 145-155. [81D.J. Filip. T.W. Ball, Procedurally representing lofted surfaces, IEEE Comput. Graphics and Applications (1989) 2733. [91E. Cohen, T. Lynche, R.F. Riesenfeld, Discrete B-spline and sub-division techniques in computer aided geometric design and computer graphics, Comput. Graphics and Image Processing 14 (1980) 137-l 11. I101W. Bohem, H. Prautzsch, The insertion algorithm, Computer Aided Geometric Design 17 (1) (1985) 58-59. Algorithm for [illE. Cohen, T. Lyche, L.L. Schumarker, degree-raising of splines, ACM Trans. Graph. 4 (3) (1985) 171-181. [I21D.F. Rogers, J.A. .4dams, Mathematical Elements for Computer Graphics, McGraw-Hill Publishing Company, New York, 1990. D31 G.V. Reklaitis, A. Ravindran, K.M. Ragsdell, Engineering Optimization - Methods and Applications, John Wiley and Sons, New York, 1984. [141J.-Y. Lai, C.-Y. Lin, Reverse engineering of composite sculp tured surfaces, Int. Journal of Advanced Manufacturing and Technology (1996) forthcoming. Cbih-Young Lin is a Ph.D. student in the Department of Mechanical Engineering at the National Central University, Taiwan, ROC. His research has focused on reverse engineering, quality control and management information system.
85 Chung-Shan Liou is a master student in the Department of Mechanical Engineering at the National Central University, Taiwan, ROC. His research has focused on CAD/CAM and reverse engineering.
Jiing-Yih Lai is a Professor in the Department of Mechanical Engineering at the National Central University, Taiwan, ROC. He obtained his M.Sc. and Ph.D. in Mechanical Engineering from the Ohio State University, Columbus, OH, USA. His research work has centered on CAD/CAM, reverse engineering and coordinate measuring metrology.