Abstract:
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.

Description:
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 1 ,K 2 ,K 3 ,K 4 , along with the power basis in the following form: 
     
         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-codes or 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&#39;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 &#34;Small errors in the approximation may result in small speed errors, it is true, but not in position errors.&#34; 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-code mean!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-code mean!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. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a diagram of the flow of information from the design system to the numerically controlled machine tool. 
     FIG. 2 shows a cubic B-spline with open end conditions. 
     FIG. 3 shows the blending functions of the cubic B-spline curve. 
     FIG. 4 shows the knot vectors for the various types of end conditions for B-spline curves. 
     FIG. 5 shows the piecewise cubic Bezier curves which represent the B-spline curve. 
     FIG. 6 shows both the B-spline control points and Bezier control points. 
     FIG. 7 shows a graphical description of the choosing of the next interpolation point on a B-spline curve. 
     FIG. 8 shows the bisection tree. 
     FIG. 9 represents the operation of bisection search. 
     FIG. 10 shows the knot vector for a B-spline curve and the parametric values and positions of the nodes of the curve. 
     FIG. 11 is a flow chart of the controlled velocity interpolation algorithm for Bezier curves. 
     FIG. 12 is a flow chart of the feed rate modification necessitated by acceleration and deceleration at the endpoints of the curve. 
    
    
     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&#39;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, γ(t), for the open or periodic control polygon B3, P={P i  } 0   n . ##EQU1## is the knot vector B2 over which B i ,k,τ (t) are defined, and B i ,k,τ (t) are the normalized local support B-spline basis functions (N1-N8 from FIG. 3) of order k(degree k-1). If τ 0  ≦τ 1  ≦. . . ≦τ 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 &#34;wrapped around&#34; 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. &#34;Unwrapping&#34; 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 &#34;Discrete B-Splines and Subdivision Techniques in Computer-Aided Geometric Design and Computer Graphics&#34; 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 &#34;A Practical Guide to Splines&#34; 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 τ=τ 0 , . . . , τ n+k , τ i  ≦τ i+1  and τ i  ≦τ i+k , for all appropriate i, and B 0 ,k,τ (t), . . . , B n ,k,τ (t) are the N+1 B-spline basis functions defined over the knot vector, τ, then for a sequence of vector or scalar coefficients P 0 , . . . , P n , comprising a control polygon and a sequence of scalars, w 0 , . . . , w n , the NURBS curve of order k is ##EQU3## 
     The {w i  } are called the homogeneous coordinates for {P 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 i  =(h x ,i,h y ,i,h z ,i, h w ,i), where if P i  =(x i ,y i ,z i ),h x ,i =w i  x i , h y ,i =w i  y i , h z ,i =w i  z i . When all the scalar values in w 0 , . . . , w 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 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 &#34;Bernstein Polynomials&#34; 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 i , (i=0, . . . , n) are the Bezier control points and θ i ,n (a, b:t) are the i th  Bernstein polynomials of degree n. The Bezier polygon is formed by joining the Bezier control points (Bezier vertices) b i , (i=0, . . . , n), in order. When the parameter t is in interval  0,1!, the Bezier curve is denoted: ##EQU8## 
     Since γ(0)=b 0 , and γ(1)=b 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 γ b 0 , . . . , b 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 &#34;A Theoretical Development for the Computer Generation and Display of Piecewise Polynomial Surfaces&#34; 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 0 ,b 1 ,b 2 ,b 3  } B6, {b 0  &#39;,b 1  &#39;,b 2  &#39;,b 3  &#39;} B7, . . . , {b 0  &#34;&#34;,b 1  &#34;&#34;,b 2  &#34;&#34;,b 3  &#34;&#34;} 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 0 ,b 1 ,b 2 , and b 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: 
     
         θ.sub.0,3 (t)=(1-t).sup.3 
    
     
         θ.sub.1,3 (t)=3(1-t).sup.2 t 
    
     
         θ.sub.2,3 (t)=3(1-t)t.sup.2 
    
     
         θ.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 ΔG4 is the parametric substep size, the initial guess of parametric substep size G4, Δ 0 .sup. 0!, for the first interpolation point is ##EQU13## where dt 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, 
     
         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 
     
         dt.sub.j.sup. i+1! =dt.sub.j.sup. i! -Δ.sub.j.sup. i+1!(9) 
    
     and operation R T3 as the iterative assignment of 
     
         dt.sub.j.sup. i+1! =dt.sub.j.sup. i! +Δ.sub.j.sup. i+1!(10) 
    
     then L T2 and R T3 operate on a segment C T1 to yield the &#34;left&#34; and &#34;right&#34; 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, α, 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 0 . 
     1. Let t j .sup. i!, t j .sup. i! =t j-1  +dt j .sup. i!, be the current guess for the parametric location of the next interpolation point. Also let Δγ j .sup. i! =∥γ(t j-1  +dt j .sup. i!)-γ(t j-1 )∥ be the chord-length distance between the last interpolation point and the current guess of the next interpolation point. When Δγ j .sup. i! is within some tolerance of the desired step size, d 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. |Δγ j .sup. i! -d 0  |&lt;ε. The parametric step size dt j .sup. i! is satisfactory. Set t j  =t j-1  +dt j .sup. i! and output γ(t j ) as the next interpolation point. Set dt j+1 .sup. 0 ! =dt j .sup. i! as the next initial parametric step size and ##EQU16## as the next initial substep size. 2. Δγ j .sup. i! -d 0  &gt;ε. The parametric step size dt j .sup. i! has caused the interpolation point to &#34;overshot&#34; the desired distance. Reduce the next parametric step size guess by the substep size, dt j .sup. i+1! =dt j .sup. i! -Δ j .sup. i!, and cut the next parametric substep size in half, ##EQU17## 3. Δγ j .sup. i! -d 0  &lt;-ε In this case, the parametric step size, dt 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 j .sup. i+1! =dt j .sup. i! +Δ 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 j+1  -t j  &lt;2dt 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 0 .sup. 0!. 
     The method of the present invention is to determine the initial parametric step size dt 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 τ={τ i  } 0   q=n+k  is a knot vector K1 for a given B-spline curve K2, then node values K3 of the curve are τ*={τ i  *} 0   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.τ,γ (t), differs from the value of the curve at the corresponding node, γ(τ i  *), by a quadratic term of t. That is, 
     
         γ(τ.sub.i *)-L.sub.τ,γ (t)=O)(|t|.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 Δs i  is the arclength step size on the curve between node points, γ(τ i  *) and γ(τ i+1  *), ΔP i  is the length of the segment of the control polygon between coefficients P i  and P i+1 , Δτ i  is the parametric step size in the interval between the node values τ i  * and τ i+1  *, and Δn i  is the parametric interval between the two nodes τ i  * and τ i+1  *. When i=0, then Δτ 0  =dt.sup. 0!, ΔP 0  =∥P 1  -P 0  ∥, Δs i  =d 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 0 , and the actual step size, d i , is determined to within some tolerance, ε, that is, 
     
         ∥d.sub.i -d.sub.0 ∥≦ε     (16) 
    
     The smaller the value of ε, 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 as 
     
         dlt=∥γ(t)-γ(1)∥t+Γ(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.τ,γ (t) and L.sub.τ,γ (t). Suppose the domain of the NURBS curve is  a, b!. If Δg i  =∥γ(τ i+1  *)-γ(τ i  *)∥ then the distance left to travel function of L.sub.τ,γ (t) can be expressed as ##EQU26## and the distance left to travel function of G.sub.τ,γ (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, 
     
         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 ε 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, α, 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 α. 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 * α 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##