Method and system for spline interpolation, and their use in CNC

An iterative spline interpolation method for a numerically-controlled machine tool device is disclosed. The non-uniform rational B-spline (NURBS) curve inputted to the numerically-controlled machine tool is interpolated with a constant step size providing the ability to maintain a controlled velocity to within a specified tolerance. In addition, the distance left to travel on the curve is obtained by a unique spline node-based approximation method providing accurate acceleration and deceleration control. The rational spline interpolation method provides significant reduction in the amount of data required to produce smoothly machined pieces while providing accurate machining of conic sections not possible by previous spline interpolation methods.

This application claims the benefit of U.S. Provisional Application No. 
60/003396, filed Aug. 30, 1995. 
BACKGROUND OF THE INVENTION 
1. Field of the Invention 
The present invention relates to numerical control of a machine tool device 
and in particular to a spline interpolation method for determining tool 
motion within a specified tolerance of the given spline curve with 
controlled velocity. 
2. Description of the Related Art 
Numerically controlled machine tools have typically used linear and 
circular arc segments along which to drive the tool. These segments are 
typically inputted to the controller by command data known as G-codes. 
Each G-code is composed of the alphabetic character G, a numerical value 
following the character G, and additional codes following the numerical 
value. More recently, part geometry has become more complex requiring the 
use of more advanced curve geometries. Parametric spline curves are 
becoming more popular in design systems and, until recently, these curves 
were approximated by a large number of linear segments before sending them 
to the controller. Unfortunately, to achieve the desired tolerance, a 
large number of very small line segments are produced. The large number of 
linear segments can be prohibitive given the memory capabilities of the 
numerically controlled machine. In addition, the shortness of the linear 
segments can cause jerky motion at high speeds, possibly even causing the 
machine to halt. To solve this problem, several spline interpolation 
methods have been developed. These method reduce the amount of data 
required by the controller and allow the controller to more adequately 
control the motion. 
Prior art relating to spline interpolation varies in the type of curves 
that can be inputted and the ability of the controller to adequately 
control the motion. Most spline functions that can be interpolated by the 
prior art are spline functions with a power basis. Such functions can be 
described by a set of coefficients, K.sub.1,K.sub.2,K.sub.3,K.sub.4, along 
with the power basis in the following form: 
EQU q(t)=K.sub.4 t.sup.3 +K.sub.3 t.sup.2 +K.sub.2 t+K.sub.1. 
Typically, a set of points to be interpolated by the spline curve is 
inputted by custom G-codesor other means!, and various conditions 
regarding smoothness and derivatives at end points are assumed or inputted 
using custom G-codes as well. The resulting piecewise spline curves are 
then evaluated at successive parametric values producing a plurality of 
positions through which the machine tool is linearly driven. 
Specifically, interpolation of a plurality of points in sequence by a 
sequence of piecewise cubic spline curves using a power basis with Hermite 
end conditions is known from U.S. Pat. No. 5,140,236 to Kawamura et al. 
Interpolation of cubic spline curves using a power basis with polynomial 
or rational Taylor series predictors is known from U.S. Pat. No. 5,321,623 
to Ensenat et al. An iterative method of spline interpolation using a 
power basis with inputted derivatives is known from U.S. Pat. No. 
4,794,540 to Gutman et al. A method for approximating an ideal path curve 
with piecewise spline curves using a power basis is known from U.S. Pat. 
No. 5,028,855 to Distler et al. All the above-cited patents are hereby 
incorporated by reference. 
However, there are three problems which are either not addressed or not 
adequately solved in the prior art 
1. Parametric curves in general, and spline curves in particular, are not 
parameterized by arclength and, except for rare special cases, such a 
parameterization can only be approximated. Without explicit control of the 
variation in the parameterization to within a given tolerance, the machine 
tool's velocity may fluctuate around the desired feed rate, possibly 
causing aberrations in the finish or unnecessary wear of the machine tool. 
2. Not all curve geometry is easily represented by splines using a power 
basis. Not only is the geometry of the power basis spline hard to 
visualize, but computer aided design (CAD) systems often do not use a 
power basis for representing spline curves. Instead, shape approximating 
splines must be converted to power basis splines before they are sent to 
the controller. In addition, splines from a power basis can only 
approximate conic sections such as circular and ellipsoidal arcs. 
3. The machine tool is very massive and must accelerate and decelerate to 
prevent gouging at the beginning and end of motion blocks in which the 
direction changes. Computing an accurate measure of how much further the 
machine tool must move along a spline curve is essential to enable the 
controller to decelerate at the appropriate time. The distance left to 
travel is non-trivial for parametric curves in general, and spline curves 
in particular, and the prior art does not address this issue. 
Maintaining constant velocity during interpolation of spline curves is 
addressed by two prior art teachings. The first, U.S. Pat. No. 4,794,540 
to Gutman et al., describes a method that requires that M-1 derivatives of 
the spline curves be inputted. The alternative to inputting these 
derivatives is to compute them by solving a system of linear equations. 
Although inputting all of the derivatives for these points is feasible for 
splines from a power basis, it is not desirable. In addition, it is not 
practical for some splines which use other basis functions, and the 
computations required to invert the matrices which are used to compute 
these derivatives are prohibitive for real-time control. Indication is 
given in Gutman that the suggested derivative approximations made to 
enforce constant velocity produce overshoot, which is unacceptable for CNC 
machining. In addition, splines which use other basis functions, such as 
non-uniform rational B-splines, are usually not derived in closed form 
because of their complexity, so the method of Gutman would not be useful 
for these types of curves. 
In the second prior art teaching, U.S. Pat. No. 5,321,623 to Ensenat et 
al., constant velocity is maintained by basing the parametric step size on 
a first or second order Taylor series approximation for the spline curve 
at each step. However, since the next parametric value is determined by an 
approximation to the derivative whereupon the position on the spline curve 
is evaluated, the actual step size may vary, causing slight fluctuations 
in the velocity of the machine tool. The Ensenat admits that "Small errors 
in the approximation may result in small speed errors, it is true, but not 
in position errors." Ensenat does not provide a method for bounding either 
the velocity errors or the deviation in actual step size, so the method 
may cause finish aberrations. 
The second problem discussed above states that the power basis is not 
sufficient to represent all curve geometry. In particular, CAD systems 
often create curves from circles and arcs which do not have exact 
representations in splines with a power basis. To represent circular 
sections, the spline curves must be rational. A particular form of 
rational spline curve which is becoming a standard type of curve used in 
computer-aided design systems is the non-uniform rational B-spline (NURBS) 
curve, described in detail later. U.S. Pat. No. 5,227,978 to Kato 
discloses a method for inputting NURBS by G-codes to a numerical control 
device, but it lacks the velocity control required for practical use in a 
numerically controlled machine tool device. In addition, the method used 
to compute the machine tool positions is recursive and, as such, is 
computationally expensive. If higher order curves are inputted, the CNC 
machine may take too long between successive outputs of position, causing 
the machine tool to stop and start. 
The third problem relates to practical embodiment of the invention. In 
order to prevent overshoot which results in gouging, the controller uses a 
measure of the distance left to travel before the end of the motion block. 
By modifying the desired step size, the velocity of the machine tool can 
not only be made constant, but can be controlled. Utilization of this type 
of acceleration and deceleration control is known for lines and arcs, the 
standard motion blocks that all controllers use, from U.S. Pat. No. 
3,727,191 to McGee. Determining the distance left to travel for parametric 
curves is more difficult but just as necessary for preventing overshoot. 
The prior art which relates to spline interpolation makes no mention of 
this practical requirement. 
SUMMARY OF THE INVENTION 
In view of the problems listed above, it is the objective of the present 
invention to provide a method and system for interpolating general spline 
curves, and in particular rational Bezier and NURBS curves of arbitrary 
order, on a numerically controlled machine tool device, using a controlled 
step size based on the desired feed rate, and providing an approximate 
measure of the distance left to travel. 
The method of this invention is to input via G-codemean!s, the requisite 
order, number of control points, knot vector, and control points which 
describe the NURBS curve. Alternately, the requisite set of control points 
for a Bezier curve is inputted via G-codemean!s. In each case, a custom 
G-code is added to the controller to enable each new type of input. The 
method of this invention uses an adaptive forward differencing technique 
to accurately produce the requested step size which, in turn, enables 
velocity control to a given tolerance. The parametric step size is 
adaptively modified from the previous parametric step size to produce an 
arclength step size to within a given tolerance of the desired arclength 
step size along the curve. 
In addition, the method of this invention also produces a estimated measure 
of the distance left to travel which enables the controller to accelerate 
and decelerate appropriately. The estimated measure of the distance left 
to travel is a smooth blend of an upper bound and lower bound of the 
actual arclength distance left to travel. 
These features of the present invention enable a numerically controlled 
machine tool device to produce a clean finish in a variety of materials 
without gouging or causing unnecessary wear on the machine tool either by 
breaking chips improperly or by making inappropriate power demands. 
These and other features of the invention will be more readily understood 
upon consideration of the attached drawings and of the following detailed 
description of those drawings and the presently preferred embodiments of 
the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
A diagram of the presently preferred embodiment of the invention is shown 
in FIG. 1. Computer I contains two modules, a CAD/CAM apparatus 2 and an 
NC post processor 3. CAD/CAM apparatus 2 provides a framework within which 
the description of the part in the form of a surface model 8 is designed. 
Surface model 8 is passed to NC post processor 3 which decomposes surface 
model 8 into motion blocks 9 comprising NURBS 10 and Bezier 11 motion 
blocks as well as the standard line and arc motion blocks and other G-code 
blocks 12. NURBS 10 and Bezier 11 motion blocks and other G-code blocks 12 
are then transferred to CNC machine tool 4. Within CNC machine tool 4, the 
first module, controller preprocessor 5, receives NURBS 10 and Bezier 11 
motion blocks and other G-code blocks 12 and transforms NURBS motion 
blocks 10, according to methods of the present invention described later, 
into Bezier motion blocks 14. In addition, other preprocessing steps not 
related to this invention are performed. Modified control blocks 13 
containing Bezier motion blocks 14, but no NURBS motion blocks 10 are 
passed to the controller 6. Controller 6 interprets modified control 
blocks 13, performing controlled velocity interpolation 16 of Bezier 
motion blocks 14 by the method of the present invention described later, 
thereby controlling the motion and operation of cutting tool 7. 
An alternate embodiment of this invention would not convert NURBS motion 
blocks 10 into Bezier motion blocks 14 in controller preprocessor 5, but 
would allow them to pass directly to controller 6 along with Bezier motion 
blocks 14 and other modified G-codes 15. Controller 6 then combines the 
steps of conversion to Bezier and controlled velocity Bezier interpolation 
16 to control the motion and operation of cutting tool 7. When controller 
6 becomes faster and can compute the required operations in the amount of 
time available between each step, this will likely become the preferred 
embodiment. 
The problems of acceleration and deceleration and controlled velocity 
discussed above are problems for any numerically controlled machine in 
which the motion of an end effector interpolates a spline. In robotics, 
for example, overshoot can cause collisions which could be damaging to the 
robot's environment or to the robot itself. In addition, velocity 
fluctuations can cause unnecessary wear on the servos which control the 
motion of the end effector. 
The presently preferred embodiment transforms descriptions of NURBS curves 
into descriptions of Bezier curves and then interpolates the Bezier curves 
with controlled velocity and other unique properties of this invention. To 
detail how the transformation takes place, a discussion of B-splines 
curves follows. 
A B-spline curve B4 shown FIG. 2 is a shape approximating curve comprising 
an order B1, a knot vector B2 and a set of vector or scalar coefficients 
which comprise a control polygon B3. Order B1 of the B-spline curve B4 and 
knot vector B2 determine a set of basis functions N1-N8 (shown FIG. 3) in 
a fashion similar to the power basis used by the prior art, which are then 
linearly combined with the coefficients of the control polygon to describe 
the shape of B-spline curve B4. In particular, the following equations 
define B-spline curve B4, .gamma.(t), for the open or periodic control 
polygon B3, P={P.sub.i }.sub.0.sup.n. 
##EQU1## 
is the knot vector B2 over which B.sub.i,k,.tau. (t) are defined, and 
B.sub.i,k,.tau. (t) are the normalized local support B-spline basis 
functions (N1-N8 from FIG. 3) of order k(degree k-1). If .tau..sub.0 
.ltoreq..tau..sub.1 .ltoreq.. . . .ltoreq..tau..sub.q is a sequence of 
real numbers, for k=1, . . . , q, and i=0, . . . , q-k, then 
##EQU2## 
If either of the denominators in Equation (3) equals zero, that term is 
defined to equal zero. 
FIG. 4 shows three common types of end conditions for the B-spline curve; 
open E1, floating E2, and periodic E3. The Bezier curve has open end 
conditions by definition. A B-spline curve has open E1 end conditions when 
its knot vector has k multiples of both the smallest and the largest 
knots. 
A B-spline has floating E2 end conditions when the k smallest and largest 
knots are simple knots, i.e. they are not duplicated. A B-spline curve is 
uniform when the knots are uniformly spaced, and non-uniform E4, 
conversely. A B-spline curve with periodic E3 end conditions has the same 
number of coefficients (control points) in its control polygon as knots, 
and has a knot vector similar to a B-spline with floating end conditions. 
In the case of floating end conditions, basis functions N1-N8 (from FIG. 
3) are non-zero outside of the parametric range of the curve. In the case 
of periodic E3 end conditions, the knot vector is "wrapped around" so that 
the curve near the beginning of the parametric range becomes a linear 
combination of control points at the beginning and end of the control 
polygon, and similarly for the end of the curve. 
Floating E2 and periodic E3 end conditions are easily converted to open E1 
end conditions. "Unwrapping" the knot vector and control polygon of a 
periodic B-spline gives an exactly equivalent floating B-spline. Adding 
knots to increase the multiplicity of both end knots E5 to k-fold (k is 
the order of the curve) converts a floating E2 B-spline into an exactly 
equivalent open B-spline. This may be done by using the Oslo Algorithm 
desribed in "Discrete B-Splines and Subdivision Techniques in 
Computer-Aided Geometric Design and Computer Graphics" by E. Cohen, T. 
Lyche and R. Riesenfeld, published in Computer Graphics and Image 
Processing, Vol. 14, 1980, which is hereby incorporated by reference. In 
the preferred embodiment of this invention, the end conditions are 
inferred from the number of control points in control polygon B3, the 
number of knots in knot vector B2, and the configuration of the knots in 
knot vector B2. 
Additional properties of B-spline curves are known, for example, from "A 
Practical Guide to Splines" by C. de Boor, published by Springer-Verlag, 
New York, 1978, which is hereby incorporated by reference. 
Further, rational B-splines have become widely used, especially to 
represent conic sections such as circular and ellipsiodal arcs. The 
non-uniform rational B-spline (NURBS) curve is described as follows. If 
.tau.=.tau..sub.0, . . . , .tau..sub.n+k, .tau..sub.i 
.ltoreq..tau..sub.i+1 and .tau..sub.i .ltoreq..tau..sub.i+k, for all 
appropriate i, and B.sub.0,k,.tau. (t), . . . , B.sub.n,k,.tau. (t) are 
the N+1 B-spline basis functions defined over the knot vector, .tau., then 
for a sequence of vector or scalar coefficients P.sub.0, . . . , P.sub.n, 
comprising a control polygon and a sequence of scalars, w.sub.0, . . . , 
w.sub.n, the NURBS curve of order k is 
##EQU3## 
The {w.sub.i } are called the homogeneous coordinates for {P.sub.i }, 
because of their resemblance in function to the homogeneous coordinates in 
a perspective transformation. The determining parameters of the rational 
B-spline curve are the knot vector and the homogeneous point h.sub.i 
=(h.sub.x,i,h.sub.y,i,h.sub.z,i, h.sub.w,i), where if P.sub.i 
=(x.sub.i,y.sub.i,z.sub.i),h.sub.x,i =w.sub.i x.sub.i, h.sub.y,i =w.sub.i 
y.sub.i, h.sub.z,i =w.sub.i z.sub.i. When all the scalar values in 
w.sub.0, . . . , w.sub.n equal 1, Equation (4) is same as Equation (1). 
The B-spline basis functions satisfy 
##EQU4## 
A Bezier curve is a polynomial curve that is a linear combination of 
Bernstein basis functions. The i.sup.th Bernstein polynomial of degree n 
over the arbitrary interval a,b! is 
##EQU5## 
The set of all Bernstein polynomials of degree n form a basis for the set 
of all polynomials of degree n. Bernstein polynomials are known, for 
example, from "Bernstein Polynomials" by G. G. Lorentz, published by 
University of Toronto Press, 1953, which is hereby incorporated by 
reference. When the interval is 0,1!, the Bernstein polynomials are: 
##EQU6## 
The representation of a Bezier polynomial of degree n is 
##EQU7## 
where b.sub.i, (i=0, . . . , n) are the Bezier control points and 
.theta..sub.i,n (a, b:t) are the i.sup.th Bernstein polynomials of degree 
n. The Bezier polygon is formed by joining the Bezier control points 
(Bezier vertices) b.sub.i, (i=0, . . . , n), in order. When the parameter 
t is in interval 0,1!, the Bezier curve is denoted: 
##EQU8## 
Since .gamma.(0)=b.sub.0, and .gamma.(1)=b.sub.n from Equation (5) the 
Bezier curve passes through first and last Bezier vertices. Also, the 
tangents at the endpoints are in the same direction as the first and last 
line segments of the control polygon. 
The Bezier curve is a special case of the B-spline curve. When knot values 
for the B-spline in Equation (3) are specially chosen, the B-spline basis 
functions reduce to the Bernstein/Bezier blending functions. Therefore, a 
B-spline curve can be exactly represented by a series of piecewise Bezier 
curves. 
Shape approximation with rational splines is generally not well understood 
except within the context of conic sections, such as arcs and ellipses. 
Indeed, rational B-splines are necessary for the exact representation of 
conic sections; ordinary B-splines or splines with a power basis are 
insufficient. 
The presently preferred embodiment incorporates a conversion from the 
general NURBS representation to the Bezier representation of the same 
order. The transformation method utilizes a known method called 
subdivision described as follows. 
Given .gamma.b.sub.0, . . . , b.sub.n :a,d!(t) as in Equation (6) then 
##EQU9## 
where k is the order and 
##EQU10## 
For a formal proof of this method of conversion, refer to "A Theoretical 
Development for the Computer Generation and Display of Piecewise 
Polynomial Surfaces" by J. Lane and R. Riesenfeld, published in IEEE 
Transactions, Vol. PAMI 2, No. 1, January 1980, which is hereby 
incorporated by reference. 
Although Bezier curves and surfaces are special cases of NURBS curves, they 
lack some of the essential features of NURBS curves. Bezier basis 
functions are global in that moving any control point will affect the 
entire curve. In addition, as the number of control points increases, so 
does the order of the curve. This can make the curve very unmanageable. 
Therefore, utilizing NURBS curves for part design will generally be the 
preferred method. 
FIG. 5 is an example of a cubic B-spline curve represented by the piecewise 
Bezier curves. B-spline curve B4 is defined by control polygon B3 and a 
uniform knot vector B2. FIG. 3 shows the blending function curves for knot 
vector B2 and order B1. FIG. 5 shows the same cubic curve shown in FIG. 2, 
described instead by its equivalent Bezier control polygons B6-B10. The 
coefficients of the control polygons {b.sub.0,b.sub.1,b.sub.2,b.sub.3 } 
B6, {b.sub.0 ',b.sub.1 ',b.sub.2 ',b.sub.3 '} B7, . . . , {b.sub.0 
"",b.sub.1 "",b.sub.2 "",b.sub.3 ""} B10 define five different cubic 
Bezier curves B11-B15. B-spline curve B4 from FIG. 2 and piecewise Bezier 
curve B11-B15 are superimposed on each other in FIG. 6 to show their 
relationship. 
In FIG. 1, when Bezier motion blocks 14 are given to controller 6, XYZ 
interpolation points 16 on the Bezier curve must be generated and given to 
machine tool device 7 at a certain update rate. The actual update rate 
will depend on the update or cycle rate of controller 6 and may be 
adjustable, but will generally be in the range of one point every 10 
milliseconds to one point every millisecond. The actual distance between 
XYZ interpolation points 16 depends on both the desired feed rate and the 
update rate with the distance being the feed rate divided by the update 
rate. To generate a smooth curve or surface, the distance between each 
interpolation point must be relatively small. To keep the distance between 
the interpolation points small with a fast feed rate, a fast update rate 
is required. Also, for machine tool 7 to run in a smooth manner with a 
constant feed rate, the spatial distance between each point must be 
relatively constant. Therefore, the problem is to generate equally spaced 
points along the Bezier curve as quickly as possible. 
For non-rational cubic Bezier curves, the curve is represented as 
##EQU11## 
The rational cubic Bezier curves is 
##EQU12## 
where b.sub.0,b.sub.1,b.sub.2, and b.sub.3 are the four vector 
coefficients of the control polygon of a cubic Bezier curve. From Equation 
(5), the basis functions of cubic Bezier curves are: 
EQU .theta..sub.0,3 (t)=(1-t).sup.3 
EQU .theta..sub.1,3 (t)=3(1-t).sup.2 t 
EQU .theta..sub.2,3 (t)=3(1-t)t.sup.2 
EQU .theta..sub.3,3 (t)=t.sup.3 
Before the first curve position is evaluated by using Equation (8), the 
initial value of the parameter t should be determined. Refer to FIG. 7 for 
the following discussion. If dt G3 is the parametric step size and 
.DELTA.G4 is the parametric substep size, the initial guess of parametric 
substep size G4, .DELTA..sub.0.sup.0!, for the first interpolation point 
is 
##EQU13## 
where dt.sub.0.sup.0! is the first guess of the parametric step size for 
the initial interpolation point discussed later. 
Since the method of the present invention is an iterative method, the last 
parametric step size is used as an initial guess of the parametric step 
size for the next interpolation point, 
EQU dt.sub.j+1.sup.0! =dt=t.sub.j -t.sub.j-1. 
For a cubic curve, the rate of change of the acceleration of the curve is 
linear. In a parameter space, this acceleration is a linear function of 
the parameter t. When the arclength step size along the curve is small, 
the change in the parameter t is small. In most cases, the last parametric 
step size is close to the desired parametric step size. If the 
interpolation point on the curve calculated by using the initial guess for 
the parametric step size is not within the tolerance of the desired step 
distance, the current parametric step size is adjusted. 
FIG. 8 is used to explain the adjustments to the parametric step size. If 
the operation L T2 is denoted as the iterative assignment of 
EQU dt.sub.j.sup.i+1! =dt.sub.j.sup.i! -.DELTA..sub.j.sup.i+1!(9) 
and operation R T3 as the iterative assignment of 
EQU dt.sub.j.sup.i+1! =dt.sub.j.sup.i! +.DELTA..sub.j.sup.i+1!(10) 
then L T2 and R T3 operate on a segment C T1 to yield the "left" and 
"right" halves, LC T4 and RC T5, of C T1 in FIG. 8. Each node in the graph 
of FIG. 8 represents a particular step size where C T1 is the initial 
guess of the parametric step size and LC T4 is half of C T1. The 
chord-length distances that the parametric step sizes LC T4 and RC T5 
represent, might not be equal. The bisection method as applied in this 
invention cuts the parametric substep in half after each iteration, 
##EQU14## 
Alternately, any factor, .alpha., can be used to scale the step size after 
each iteration. The formula would then be 
##EQU15## 
The iteration terminates when the interpolation point evaluated at some 
node in FIG. 8 is within the specified tolerance of the desired step 
distance, d.sub.0. 
1. Let t.sub.j.sup.i!, t.sub.j.sup.i! =t.sub.j-1 +dt.sub.j.sup.i!, be 
the current guess for the parametric location of the next interpolation 
point. Also let .DELTA..gamma..sub.j.sup.i! =.parallel..gamma.(t.sub.j-1 
+dt.sub.j.sup.i!)-.gamma.(t.sub.j-1).parallel. be the chord-length 
distance between the last interpolation point and the current guess of the 
next interpolation point. When .DELTA..gamma..sub.j.sup.i! is within some 
tolerance of the desired step size, d.sub.0, the search is stopped and the 
forward operation is made. There are three possible cases that need to be 
considered for each successive guess until a satisfactory guess is 
reached. 
1. .vertline..DELTA..gamma..sub.j.sup.i! -d.sub.0 .vertline.&lt;.epsilon.. 
The parametric step size dt.sub.j.sup.i! is satisfactory. Set t.sub.j 
=t.sub.j-1 +dt.sub.j.sup.i! and output .gamma.(t.sub.j) as the next 
interpolation point. Set dt.sub.j+1.sup.0 ! =dt.sub.j.sup.i! as the next 
initial parametric step size and 
##EQU16## 
as the next initial substep size. 2. .DELTA..gamma..sub.j.sup.i! -d.sub.0 
&gt;.epsilon.. The parametric step size dt.sub.j.sup.i! has caused the 
interpolation point to "overshot" the desired distance. Reduce the next 
parametric step size guess by the substep size, dt.sub.j.sup.i+1! 
=dt.sub.j.sup.i! -.DELTA..sub.j.sup.i!, and cut the next parametric 
substep size in half, 
##EQU17## 
3. .DELTA..gamma..sub.j.sup.i! -d.sub.0 &lt;-.epsilon. In this case, the 
parametric step size, dt.sub.j.sup.i!, is not large enough to cause the 
next interpolation point to be within range of the desired distance. 
Increase the next parametric step size guess by the substep size, 
dt.sub.j.sup.i+1! =dt.sub.j.sup.i! +.DELTA..sub.j.sup.i!, and cut the 
next parametric substep size in half, 
##EQU18## 
This particular embodiment assumes that there is a limit on the number of 
iteration steps that can be taken to compute a given interpolation point. 
It also assumes that the next parametric step size will be less than twice 
the size of the current step size, i.e. t.sub.j+1 -t.sub.j &lt;2dt.sub.j. If 
the interpolation point at twice the current step size is not far enough 
along the curve to be within the tolerance of the desired step size, the 
next parametric step size halves the distance to twice the current step 
size until the number of iterations exceeds the limit. 
An example is described by FIG. 9. Initial step size L1 is added to the 
parametric position of the last interpolation point L2 to obtain the 
initial guess of the parametric location of the next interpolation point 
L3. After computing the chord-length distance between the last 
interpolation point and the initial guess L3, it is determined that 
initial guess L3 overshoots desired step size L13. Case 2 above is 
utilized and current substep size L4 is subtracted from parametric step 
size L1, to produce second parametric step size L5. Second parametric step 
size L5 is added to the parametric position of last interpolation point L2 
to obtain second guess L6 of the parametric location of the next 
interpolation point. After the chord-length distance to second guess L6 is 
computed, it is determined that second guess L6 undershoots desired step 
distance L13. Case 3 from above then applies, and second substep L7 is 
added to second parametric step size L5 to produce third parametric step 
size L8. Once again, after the chord-length distance to third guess L9 is 
computed, it is determined that third guess L9 undershoots desired step 
distance L13. Case 3 from above applies again and third substep L10 is 
added to third parametric step size L8 to produce fourth parametric step 
size L11. When added to the parametric location of the last interpolation 
point L2, fourth guess L12 is within accepted tolerance L14 of desired 
step distance L13. Case 1 applies from above and fourth guess L12 is 
outputted to machine tool 7. 
For each cubic Bezier segment of the curve (either rational or 
non-rational), the parameter t is defined from 0 to 1. Even though the 
bisection method as applied above can adjust t up or down very quickly, it 
is still necessary to find a good initial value for dt for the first 
interpolation point to minimize the adjustment time. Given an initial 
guess for dt in (0, 1), the adaptive method may adjust down or adjust up 
many times until the step is found to be within the tolerance to do the 
forward operation. If the calculation is done outside CNC controller, the 
initial scaling down of dt may not be required. If, however, the 
calculation is done the inside CNC controller, the desired step L13 must 
be found within the update rate. It is therefore necessary to determine a 
reasonable initial parametric step size, dt.sub.0.sup.0!. 
The method of the present invention is to determine the initial parametric 
step size dt.sub.0.sup.0! based on properties of B-splines relating to 
node values of the B-spline curve. Refer to FIG. 10 for the following 
discussion. If .tau.={.tau..sub.i }.sub.0.sup.q=n+k is a knot vector K1 
for a given B-spline curve K2, then node values K3 of the curve are 
.tau.*={.tau..sub.i *}.sub.0.sup.q=n where 
##EQU19## 
Since control polygon K4 has a direct correspondence to nodes N2, control 
polygon K4 can be made into a piecewise linear curve with the parametric 
values at the coefficients of the control polygon being equal to 
corresponding nodes K5-K10 according to the following equation: 
##EQU20## 
It is known, then, by the variation diminishing property of B-splines, 
that the linearized polygon, L.sub..tau.,.gamma. (t), differs from the 
value of the curve at the corresponding node, .gamma.(.tau..sub.i *), by a 
quadratic term of t. That is, 
EQU .gamma.(.tau..sub.i *)-L.sub..tau.,.gamma. 
(t)=O)(.vertline.t.vertline..sup.2) (13) 
Therefore, linearized polygon K4 is a good approximation to curve K3 at 
node points K5-K10. From this, a relationship between the arclength of the 
curve and the parametric domain of the curve can be established. In 
particular, 
##EQU21## 
over a given interval, i, where .DELTA.s.sup.i is the arclength step size 
on the curve between node points, .gamma.(.tau..sub.i *) and 
.gamma.(.tau..sub.i+1 *), .DELTA.P.sup.i is the length of the segment of 
the control polygon between coefficients P.sub.i and P.sub.i+1, 
.DELTA..tau..sup.i is the parametric step size in the interval between the 
node values .tau..sub.i * and .tau..sub.i+1 *, and .DELTA.n.sup.i is the 
parametric interval between the two nodes .tau..sub.i * and .tau..sub.i+1 
*. When i=0, then .DELTA..tau..sup.0 =dt.sup.0!, .DELTA.P.sup.0 
=.parallel.P.sub.1 -P.sub.0 .parallel., .DELTA.s.sup.i =d.sub.0, and from 
Equation (11), 
##EQU22## 
for the Bezier case. Solving for dt.sup.0! by substitution into Equation 
(14), 
##EQU23## 
The method of the present invention utilizes an adaptive forward difference 
operation to determine the next step size, since the relationship between 
the arclength of the curve and the parameterization of that curve is not 
constant. The difference between the desired arclength step size, d.sub.0, 
and the actual step size, d.sub.i, is determined to within some tolerance, 
.epsilon., that is, 
EQU .parallel.d.sub.i -d.sub.0 .parallel..ltoreq..epsilon. (16) 
The smaller the value of .epsilon., the more adaptive steps are taken to 
determine the next step size. 
The desired arclength step size is inputted to the adaptive method 
described here and is utilized to control the velocity at which the 
machine tool device travels. To determine the desired arclength step size, 
controller 6 requests the distance left to travel from the interpolation 
module so that it can accelerate and decelerate appropriately. A simple 
embodiment of this method of the present invention takes smooth blend of 
the distance between the current position on the curve and the end of the 
curve and the distance left to travel computed from an arclength 
parameterization of the control polygon scaled to match the parametric 
domain of the curve. For a Bezier curve, the parametric domain is 0,1!. 
Suppose that 
##EQU24## 
is the normalized distance left to travel function for the control 
polygon. Then the simple embodiment computes the distance left to travel a 
s 
EQU dlt=.parallel..gamma.(t)-.gamma.(1).parallel.t+.GAMMA.(t)(1-t).(17) 
Although this is clearly an approximation, as the machine tool approaches 
the start of deceleration, the approximation becomes more accurate and 
produces acceptable results. This method works best in the context of 
Bezier curves since the parameterization is well-behaved. 
An alternate method for computing the distance left to travel that produces 
more accurate results for NURBS curves in particular, is the following. 
Along with the piecewise linear function defined by Equation (12), another 
piecewise linear function is defined by connecting the node points of the 
curve with a similar parameterization to that of Equation (12). Thus, 
##EQU25## 
In this method the distance left to travel is determined by an average of 
the distance functions of the two curves, G.sub..tau.,.gamma. (t) and 
L.sub..tau.,.gamma. (t). Suppose the domain of the NURBS curve is a, b!. 
If .DELTA.g.sup.i =.parallel..gamma.(.tau..sub.i+1 *)-.gamma.(.tau..sub.i 
*).parallel. then the distance left to travel function of 
L.sub..tau.,.gamma. (t) can be expressed as 
##EQU26## 
and the distance left to travel function of G.sub..tau.,.gamma. (t) can be 
expressed as 
##EQU27## 
The computation of dlt from these functions is 
##EQU28## 
This distance left to travel function is not particularly time consuming 
to compute and provides a substantially more accurate measure of the 
distance left to travel for NURBS curve. 
The diagram of the method of the present invention for controlled velocity 
interpolation of Bezier and NURBS curves is presented in FIG. 11 and FIG. 
12. This method corresponds to controller 6 output of interpolated 
positions 16 to machine tool 7 in FIG. 1. For the discussion of the 
following steps of this adaptive method, refer to FIG. 11 and FIG. 12. The 
input to controller 6 comprises the desired Feed Rate, the lower bound of 
the processing time per loop step, or the Loop Time, and the second 
through fourth coefficients of the control polygon of the Bezier curve to 
be interpolated. The last position of the machine tool, LP, is provided by 
the current state of controller 6. The Feed Rate and the Loop Time are 
used to compute the desired distance per step, dist, in step S1, according 
to the following formula, 
EQU dist=Feed Rate * Loop Time 
where Loop Time is proportional to the clock rate of processor 6 in the 
preferred embodiment of the invention. The Loop Time can be determined 
empirically, or can be measured by the average number of instructions 
executed per step within the loop. 
In step S2, the initial value for step is determined by the node method of 
Equation (15) using the inputted coefficients, LP, and dist. The parameter 
t is initialized to the value of step. In addition, Accelerating is set to 
TRUE and Decelerating is set to FALSE. In step S3, the next interpolated 
position, NP, is computed from Equation (8). The next position, NP is then 
evaluated to see if it is within .epsilon. of the desired distance to 
travel, dist, or if the number of iterations exceeds the limit, in step 
S4. If neither of these conditions are true, we enter step S11. Otherwise 
we proceed to step S5. 
In step S11, the difference between NP and LP is checked to see if it 
exceeds dist. If so, step is subtracted from t in step S12. Otherwise, 
step is added to t in step S13. Then in step S14, step is divided by 2 or 
another factor, .alpha., and we proceed back to the top of the loop at 
step S3. In an alternate embodiment of the present invention, the next 
value of t can be can be determined by one of several techniques. In the 
first of the alternate methods, the amount of overshoot or undershoot is 
determined by a ratio of the difference between NP and LP and the desired 
step distance, dist. This ratio is then used to determine the next value 
of t. The step size is set to be the difference between the next value of 
t and the last value of t, and then is divided by the factor of 2 or 
.alpha.. The second alternative method fits a quadratic curve to the 
distances between the next position NP, and the last two computed 
positions and determines the next parametric value of t by using this 
curve as an quadratic estimate of the distance. 
When step S5 is reached, NP is output to controller 6. In step S6, we 
proceed to step A1 of FIG. 12 to modify dist for acceleration and 
deceleration if necessary. Upon return from the process detailed in FIG. 
12 at step A10, we proceed to step S7. In step S7, we obtain a new 
parametric step size, step, by taking a ratio of the new curve length step 
size, dist to the old curve length step size, and multiplying it by the 
old parametric step size. That is, 
##EQU29## 
In step S8, t is incremented by step * 2 or step * .alpha. and LP is set to 
the value of NP. Step S9 determines if the interpolation is finished by 
checking if the parametric value of t exceeds the bounds of the Bezier 
curve which is 1. If not, we proceed to step S3. If so, the interpolation 
terminates at step S10. 
FIG. 12 details the process of controller 6 which determines if the 
velocity needs to increase or decrease according to the parametric 
position t on the inputted curve. Step A1 of FIG. 12 is entered from step 
S6 of FIG. 11. Step A2 computes the distance to stop, dts, at the current 
Feed Rate, Acceleration and desired step size, dist. Also in step A2, the 
distance left to travel, dlt, is computed from Equation (17) or Equation 
(21). In step A3, the controller determines if it should be accelerating 
or decelerating by setting Decelerating to true if the distance left to 
travel, dlt, minus the desired step size, dist, is less than the distance 
to stop, dts, and by setting Accelerating to TRUE if Decelerating is FALSE 
and the desired step size is less than the largest step size, max, allowed 
by the controller 6. Step A4 checks if controller 6 should decelerate. If 
so, controller 6 proceeds to step A5. If not, controller 6 skips to step 
A7. 
In step A5, dist is decremented by the Acceleration provided by the 
controller 6. In step A6, if dist is smaller than the smallest step size 
allowed by controller 6, min, then dist is set to min. Controller 6 then 
proceeds to step A7. 
In step A7, if controller 6 determines that it is not Accelerating, then it 
returns to step S6 of FIG. 11. Otherwise, dist is incremented by 
Acceleration in step AS, checked against the largest step size allowed by 
controller 6, max. If dist is greater than max, in step A9, dist is set to 
max in step A12. In step A10, controller 6 returns to step S6 of FIG. 11. 
One realized embodiment of the invention for testing CPU time was written 
in C language running in UNIX with a 68030 microprocessor. Computer source 
code which implements this embodiment is in Appendix A. The results shown 
that the average CPU run time is about 150 microseconds for each loop for 
the adaptive method and about 240 microseconds for the subdivision method. 
Therefore, the approximate maximum CPU time required to calculate the next 
position on the cubic Bezier curves with an equal forward step is about 1 
millisecond using the method of the present invention (running on a 68030 
processor). 
Computer programs implementing the method of this invention will commonly 
be distributed to users on a computer-readable medium such as floppy disk 
or CD-ROM. From there, they will often be copied to a hard disk or a 
similar intermediate storage medium. When the programs are to be run, they 
will be loaded either from their distribution medium or their intermediate 
storage medium into the execution memory of the computer, configuring the 
computer to act as a mass storage emulator or mass storage access program. 
All these operations are well-known to those skilled in the art of 
computer systems. 
It is to be understood that the above described embodiments are merely 
illustrative of numerous and varied other embodiments which may constitute 
applications of the principles of the invention. Such other embodiments 
may be readily devised by those skilled in the art without departing from 
the spirit or scope of this invention and it is our intent they be deemed 
within the scope of our invention. 
##SPC1##