ERROR-FREE INTERPOLATION OF PARAMETRIC SURFACES J.G. Gan Nanyang Technological Institute Singapore T.C. Woo Dept. of Industrial & Operations Eng. The University of Michigan Ann Arbor, MI 48109-2117 Technical Report 89-35 December 1989 Revised February 1991

ERROR-FREE INTERPOLATION OF PARAMETRIC SURFACES* J.G. Gan T.C. Woot Nanyang Technological Institute University of Michigan Singapore 1205 Beal Ave. Ann Arbor, MI 48109-2117 December 12, 1989 Revised: February 1991 Abstract A Digital Differential analyzer (DDA) is a cost effective implementation for interpolating parametric curves and surfaces. However, the problems of register overflow and integer round-off have prevented it from being adopted by industry. Contrary to intuition, an increase in register capacity is shown to have no bearing on the register overflow problem. It is shown that there is a minimum number of interpolation steps Ke below which register overflow occurs. Correspondingly, it is also shown that there exists a maximum number of interpolation steps KU above which the accumulation of integer round-off errors exceeds the acceptable limit. But, when Ke > K,, conflict arises. A solution is given to resolve the possible conflict and to yield an error-free interpolation of curves and surfaces. *To appear in ASME Transactions, J. of Engineering for Industry tPlease send correspondence to T.C. Woo at the above address 1

1 INTRODUCTION A numerical control (NC) machine executes according to a part program which is converted into a set of low-level instructions. One of the instructions is to pulse the axis motors such that the cutting tool moves, with respect to the workpiece, along a specified path by interpolation. The available NC interpoators have been linear and circular. Under linear interpolation, up to three linear axes move simultaneously such that a straight line results [Milner 76, Koren 76]. The circular interpolator is a little more restrictive. Only two linear axis motors are pulsed such that the resulting arc lies in one of the three principal planes [Koren 81]. More recently, some of the NC machine controller manufactures are beginning to offer spline interpolation, by first approximating a spline with a series of line segments, as illustrated in Figure la. The error e of such an approximation is a function of the maximum length of the line segments e and the minimum radius of curvature p [Faux 79]: e 8p Since the resolution of modern NC machines has been improved to as high as 0.1 gm, it is reasonable to expect shorter line segment lengths for spline interpolation. Indeed, the NC machine controller manufacturers, with the aid of high-speed 32-bit microprocessors and multi-master-bus architecture, offer linear interpolation on the linear approximation of a spline, as shown in Figure lb. This paper offers error-free interpolation of spline curves and surfaces, as illustrated by Figure Ic. We say that the interpolation is error-free when the error is no greater than 2 of the available resolution of the NC machine. (a) Approximation (b) Interpolation (c) Error-free of Approximation Interpolation Figure 1 Approximation and Interpolation 2

It may be asked: Why don't the NC manufacturers approximate a curve to the resolution of the machine, hence achieving error-free interpolation by effectively by-passing linear approximation? Clearly, the error is dominated by approximation. As attempts of non-success are seldom reported, one can only surmise that there are difficulties. In this paper, we report the difficulties that have prevented successful implementation of error-free interpolation of curves and surfaces: overflow and round-off errors, as illustrated by Figure 2. We further present a solution that overcomes these two kinds of errors, using the Digital Differential Analyser as implementation. Register Round-off Error-free Overflow Errors Interpolation Figure 2 Interpolated Surface The Digital Differential Analyzer (DDA) scheme can be used to perform interpolation of integral, differential, exponential, trigonometric and polynomial functions, etc. Previous work [Mayorov 64, Sizer 68] detailed how the interpolation of these functions could be implemented using hardware. Researchers [Koren 76, Milner 76, Koren 81] have introduced software implementation of the DDA technique for NC machine controller applications such as line aAd arc interpolations. In this paper we examine the use of DDA for interpolating 3

curves and surfaces of arbitrary degrees. We assume that a curve is represented parametrically as a polynomial in u. X(u) = ao + alu + a2u2 +.. + aui... +. amum Y(u) = bo+bl+u2 +...+ biui +...+bmUm Z(u) = Co + clu + C22 +.. + Ci+... + CmU Stepping along the parameter u E [0,1] traces the curve. A surface of two parameters, u and v, involves keeping one of the parameters constant while stepping along the other. X(u,v) = aoo + alou +...+ aioui +...+ aolv +... + aiju'vj +.. + ammumvm Y(u,v) = boo + bou+... + bou +...+ bov+... + bijuiv+ +... +bmmumvm Z(u, V) = C+ C10U+... +CiU' +... + C01V + -.. + CijU'V3 +.. + CmmUmVm Rational curves and surfaces are represented and traced similarly. Pushing the technology of DDA beyond lines and arcs offers the potential of high speed curve and surface interpolation with inexpensive hardware for numerical controllers and for computer graphics workstations. However, this is not without its difficulties. The basic principle behind the DDA parametric interpolation scheme is the approximation of integration by additions. It can be viewed as a discrete version of the forwarddifferencing method [Ding 87, Foley 82], implemented using finite-capacity registers for data storage. Each interpolation step of a polynomial of degree m consists of the following successive addition operations: FOR i: = 1 TO m DO (Rmi) - (Rm-_) + (Rm-i,+) (1) where (Rm_i) represents the content of register Rm_i. Figure 3 provides a schematic of the implementation of (1) for interpolating a cubic curve (m = 3) using four registers R3, R2, R1 and Ro. The output of the interpolation for a polynomial function is in the form of overflow of the last register, Ro, in the cascade. It is used directly to instruct the motor to move a step in an open-loop control system, or to signal an error in the desired position in a closed-loop control system. Real-time computational speed achievable by DDA to perform the steps outlined in (1) is attributed to the usage of fixed-point variables and simple arithmetic operations such as increment, copy, compare and detect-overflow. 4

An examination of (1) reveals two problems in the implementation of the DDA parametric interpolator. Firstly, no register, other than Ro should be allowed to overflow because if any of them does, the value stored is no longer valid. (An overflow results when the value stored exceeds the capacity.) Secondly, errors due to integer round-off of register values accumulate and propagate at an exponential rate [Kanatani 84]. These two problems were illustrated in Figure 2. This paper is organized as follows. Section 2 discusses the problem of register overflow and recommends a solution to prevent it. Section 3 discusses the problem of integer round-off error and how the error can be controlled to within tolerance. Section 4 provides an algorithm to satisfy simultaneously the conditions for preventing register overflow and limiting roundoff error. We conclude in Section 5. Copy & Add nr Copy & Add Overflow Signal t - I Figure 3 R2 R1 RO Schematic Drawing of DDA 5

2 REGISTER OVERFLOW The conditions under which the registers overflow is first discussed. While intuition suggests that overflow is related to the capacity of a register, it is shown, however, that the register overflow problem cannot be overcome by increasing the capacities of registers. Specifically, there exists a lower bound on the number of interpolation steps which prevents register overflow. 2.1 Conditions for Overflow The output from the DDA parmetric interpolator is the overflow signal from the last register, Ro, in the cascade. None of the other registers is allowed to overflow, otherwise the value stored would be incorrect, resulting in wrong interpolation results. A DDA interpolates by performing successive copying of the content of register Ri, denoted by (Ri), and adding it to (R,-1), the content of register Ri-1, as noted from (1). (R,_,) is always greater than (Ri) for all i > 1. Hence, it follows that the register most likely to overflow is R1. To prevent it from overflowing, we have to ensure that (Rl)max < 2n where n is the number of bits in R1. Converely, the capacity of R1 must be n > log2 (R)max. 2.2 Register Capacity It appears that a large n will prevent R1 from overflowing. But, the following lemma shows that whether the R1 register will overflow is independent of n. Lemma 2.1 Whether register R1 overflows is independent of its capacity n. Proof: Let cra be the initial value for register Ri. The content of R1 after j integration steps is given by j(j + 1) j(j+ 1)(j+ 2) (j + 1)(j + 2)(j + 3) 2! 3 4+ h(Re values of te2+ eir divatis gin + + Append (2)I The values of the a's, their derivations given in Appendix I, are: 6

5!2" = K5D (as +D) 4!2n a4 = ^-ir^044 -- 4 =K4 * D (a4 +. ) 3!2n = K3*D(a3+' ) 2!2" 3a3 a2 =(a2- +T ~ K2 * D K 2n a2 a3 al =(a + K *D K(a K2+K.. where D = unit distance moved due to a single output pulse from DDA, and K = total number of steps to interpolate a curve. Substituting the values of the a's into (2), we have 2n a2 a3.2 a4 as (R1) = K. D{a+ (2j-1)+K2(32-3j+1)+K (..)+ K4(...).+.} = 2 * f(ai,j,K, D) (3) where 1 a2 a3.2 a4 f(ai,j,K, D) = K {(al + (2j - 1) + (3j2 - 3j + 1) + (...) + 4 (...) +..}. (4) R1 overflows when (R1) > 2n. From (3), it follows that in order to prevent R1 from overflowing, f(a, K, D) < 1, or 1 a2 as3 2 a4 as K ~ D {al + f(2j - 1) + V 2(3j2-3 + 1) + (..) + (..) + (5) K 1 D (5) < Ifor all j E [O,K]. Since f(.) is independent of n, whether R1 register will overflow is independent of n. This completes the proof. * 2.3 Analysis of Overflow The proof of Lemma 2.1 shows that the overflowing of a register is a function of K, the number of steps for interpolating a curve the coefficients of the curve ai, as well as D, the resolution of the machine tool. D and ai are fixed. It is, therefore, interesting to investigate the lower bound on K so as to avoid the register overflow problem. Determination of such a lower bound for interpolating a polynomial funciton of any degree m is in order. 7

Lemma 2.2 There exists a minimum number of interpolation steps IKe above which register overflow does not occur. Proof: It suffices to consider the condition to prevent register R1 from overflowing as given by (5), reproduced below. K D al + (2j -1)+ 32(3j2 3Ji 1)+ (...)+.. < 1 for all j E [0, K]. (5) From (5), we seek a relation between K and the given variables D and aj. Since each of the terms can be bounded as follows: a(2j - 1) < 2a2 if a2 > 0 > 2a2 otherwise 2 (32- 3j + ) < 3a3 if a3 > O > 3a3 otherwise and so on for all j E [0, K], we have POS(al) + 2 POS(a2) +... + m POS(am) < K * D NEG(ai) + 2 NEG(a2) +... + m NEG(am) > -K * D where POS(x) = x if x > 0 = 0 otherwise and NEG(x) = -x if x < 0 = 0 otherwise A tight condition for bounding K is 1 K > - MAX {[POS(ai) + 2POS(a2) +... + mPOS(am)], D -[NEG(al) + 2NEG(a2) +... + mNEG(am)]} (6) Let K, denote the term on the right of (6), we have K = MAX { [ iPOS(ai)], iNEG(a,)] (7) which, is bounded. This completes the proof that there exists a minimum number of steps K, above which no register overflow occurs. * 8

3 INTEGER ROUND-OFF The source and propagation of round-off errors are now discussed. With the aid of numerical experiments, the control of such errors to within the tolerance limits is then established in Section 3.2. 3.1 Source and Propagation of Errors Since the registers have finite capacities, there are representation errors. As integers are stored and the fractional parts are lost, there are the round-off errors at initialization. During each interpolation step, these round-off errors accumulate and propagate. Let ci be the integer round-off error in register R, when the register is initialized. After K steps, the error in Ro register is given by K(K + 1) K(K+1)(K + 2) K(K+1)...(K+m-1) (8) eo + Kq +.. 2 2 + (3 +.. + m ( 2! 3! m, where m is the degree of the curve. This well-known problem of exponential growth of cumulative error with the number of incremental steps illustrates the major weakness of conventional incremental approach [Kanatani 84]. Since every 2n in Ro corresponds to a single output pulse, the number of pulses ~ erroneously sent out due to the initial register round-off errors is given by + K(K + 1) K(K + 1)(K + 2) = [o +K+ K2i +'6 +....+ 2! 3! K(K + 1)...(K + m- 1) Em]*2-" (9) We discuss in the next section the control of C. 3.2 Controlling Round-off Errors In Appendix I, the derivation of the initial value, a,, to be stored in register Ri shows that ai is a function of the coefficient of the curve a's, number of interpolation steps K, and the register size n. Since e, is the error due to the integer round-off of ca;, c, is also a function of a's, K and n. From (9), we deduce that C is also a function of a's, as well as K and n. 9

I10 I I I I I I>K -100 —- -L -- -I -- I - K o10 1 10 1083 104 10' 10 1 Figure 4 Round-off Errors Consider a specific cubic function x(u) = 00u + 100u2 + 100u3 (10) Figure 4 shows the number of erroneous pulses C verses the number of interpolation steps K for registers of sizes 8-bit, 16-bit, 32-bit, 64-bit and 128-bit. The plots suggest that there is an upper limit to the number of steps that can be used for DDA interpolator with different register capacities. In the case of interpolating the function (10), the upper limits for K are approximately 14, 75, 4073, 6760830 and some value greater than 109 when the sizes of the registers are respectively 8-bit, 16-bit, 32-bit, 64-bit and 128-bit. We seek a theoretical upper limit for K. Lemma 3.1 There exists a maximum number of interpolation steps Ku above which the accumulation of integer round-off errors exceeds the acceptable limit. Proof: Suppose 2 is the acceptable number of erroneous pulses ~. Since e is due to integer round-off, it may assume any value between 0 and 2. In the worst case, every e is 2. By substituting 2 into every E, and bounding C to -, the condition given by (9) becomes 1+ K K(K + 1) K(K 1)(K 2) K(K+ 1)...(K + m-1) 1 2! + +" -- r + -+ --- <. 1 2! 3! M.,I 10

Dividing both sides by K gives K+1 K2+3K+2 Km"- + rm(m+1)Km-2 +.. + (m-1) (2" - 1) 1+ + +...+ 2 2! 3! m! K Rearranging in decreasing order of K yields K"-' m(m +) 1 1 1 (2n - 1)m, -+ K-2+... + (...)K + (1 +- +...+ + -)m < or 2 2 3 m, Km-i. m(m + 1) (...) I + 1 1 1 m! (2n - l)m! K [1+ ^I^~+..(+'" +(1+ + + —-+-)m x]< -' (11) 2K Km-2 2 3 m Km-i] Since K is usually very large (in the order of 105) compared to m (the degree of the curve), we can safely assume that the first term in (11) is much greater that the second term, and so on. The terms within the square brackets in (11) sum to less than 2. The condition to restrict round-off error to 0.5 of an output pulse can therefore be simplified to (2K - 1)m! Km < 2 Therefore, the maximum number of integration steps Ku is given by taking the mth root, K=( ((2n 2 1)m! )/m (12) Since n, the number of bits, and m, the degree of the curve, are bounded, Ku is finite. In order to limit the number of erroneous pulses, care must be taken to ensure that the number of interpolation steps used is lower than Ku. If this upper limit is too low, overflow (as discussed in Section 2) can still occur. This is the subject of the next section. 11

4 COMBINED SOLUTION Lemma 2.2 asserts that, to prevent overflow, the number of interpolation steps K must exceed a lower bound 1 (Fm ] r m l 1 K>Kt, MAX iPOS(ai)], -E iNEG(ai) (7) DI Li= J L <i=1 and Lemma 3.1 states that, to prevent round-off errors, K must not exceed an upper bound K < K = (2 - l)m! )/ (12) For error-free interpolation, both conditions must be satisfied simultaneously, i.e., Ke < K < Ku: 1 {[m ] [ m ]} )(V2n )M!1/ -MAX iPOS(ai) \ iNEG^i)\ < 2 —Li L=l J L i=1 / The necessary number of bits n to ensure no overflow and no more than 4 = 2 of an erroneous pulse is given by n > -logm! + mlog [ MAX { [ iPOS(ai)], - iNEG(ai)]}] (13) since 2" > 1. If ( is., then the "1" on the right hand side of (13) is replaced by q. With condition (13), we are ready to summarize our findings in a procedural form. PROCEDURE DDA - CURVE Step 1 {Determine the minimum K that ensures no register overflow problem} m m m K - *MAX { iPOS(ai,), iPOS(ai,), iPOS(ai,), i=1 i=1 i=1 1m1 n m - iNEG(ai,),- iNEG(a,), - E iNEG(a(,)} iml i=l i=l Step 2 {Find n} n - 1 - logm! +m log K Step 3 {Calculate initial register values} FOR each axis DO (Ro) - ROUND (2n/D * (ao - D * TRUNC (ao/D))) (R1) - ROUND (2n/(K * D) * (a - a2/K + a3/K2)) 12

(R2) ROUND (2! 2n/(K2 * D) * (a2- 3a3/K)) (R3) - 3! 2n/(K3 *D) * a3 ENDDO Step 4 {Position motor} Move motor to position (D * TRUNC (ao/D)) Step 5 {Generating Output Pulses} For i:= TO K DO Wait for a system clock pulse FOR each axis DO (R2) -(R2)+ (R3) (R1) - (R) + (R2) (Ro) (Ro) + (R1) IF Ro overflows THEN Send an output pulse to motor ENDDO ENDDO ENDDO {ELSE} ENDPROCEDURE PROCEDURE DDA - Surface Approximate surface by a set of constant parametric curves For each curve DO Call DDA-Curve ENDPROCEDURE It may be noted that Procedure DDA-curve has two phases: a preprocessing phase (involving Steps 1 to 3 that need to be done only once for a given curve) and the actuation phase (Step 5 which is repetitive and arithmetically simple). The preprocessing phase can be done in software and the actuation phase implemented in hardware. The interpolated surface in Figure 2 was generated entirely by software. The bicubic equations used for generating the surfaces plotted are as follows: x(u,w) = 100 + 50 u + 100 u2 + 100 u3 + 0 w + 50uw ++ 5 u2w - 100 U3w + 50 2 - 50 uw2 + 0 u2w2 - 200 u3w2 - 50 W3 + 50 uw3 + 0 u2W3 + 100 U3w3 13

y(u,w) = 50 + 100 u + 50 u2 + 50 u3 + +0 + 50 u 0 u2w + 200 u3w + 200 w2 - 50 uw2 + 0 u2w2 - 100 u3w2 - 50 w3 + 100 uw3 + 50 u2w3 - 100 u3w3 z(u,w) = 100 + 50 u + 100 u2 + 100 u3 + 0 w + 50 uw + 50 u2w - 100 u3w + 50 w2 - 50 uw2 + 0 u2w2 - 200 u3w2 - 50 w3 + 50 uw3 + 0 u2W3 + 100 u3w3 5 SUMMARY The DDA has been shown to be an efficient scheme for interpolating not only lines [Milner 76] and arcs [Koren 81], but also parametric curves and surfaces of arbitrary complexity. However, because of the complexity, the two sources of errors, overflow and round-off, can no longer be overlooked. It is shown in this paper that overflow has nothing to do with the register size n. It is also shown the round off error is a function of n (as well as m, the degree of the curve or surface). To overcome both sources of error, the choice of n is rationalized as the simultaneous satisfaction of two constraints on the number of interpolation steps K. 14

REFERENCES 1. Ding, Q. and B.J. Davis, "Surface Engineering Geometry for Computer-Aided Design and Manufacture," Ellis Horwood, UK, 1987. 2. Faux, I.D. and M. Pratt, "Computational Geometry for Design and Manufacture," Halsted Press, NY, 1979. 3. Foley, J.D., and A. VanDam, "Fundamentals of Interactive Computer Graphics," Addison-Wesley, Reading, MA, 1982. 4. Kanatani, K., "Errors of the Incremental Method for Curves," Comp Vis Grap I Porc, Vol. 26, pp. 130-133, 1984. 5. Koren, Y., "Interpolator for a Computer Numerical Control System," IEEE Trans. on Computers, Vol. C-25, No. 1, pp. 32-37, Jan. 1976. 6. Koren, Y., and 0. Masory, "Reference-Pulses Circular Interpolators for CNC Systems," Trans. ASME, J. Eng. Ind., Vol. 103, No. 1, pp. 131-136, Feb. 1981. 7. Mayorov, F.V., "Digital Differential Analyzers," liffe Books, London, England, 1964. 8. Milner, D.A., "Some Aspects of Computer Numerical Control with Reference to Interpolation," Trans. ASME, J. Eng. Ind., pp. 883-889, Aug. 1976. 9. Press, W.H., B.P. Falnnery, S.A. Teukolsky and W.T. Vetterling, "Numerical RecipesThe Art of Scientific Computing," Cambridge U. Press, Cambridge, England, 1986. 10. Sizer, T.R., "The Digital Differential Analyser," Chapman & Hall, London, England, 1968. 15

APPENDIX I Calculation of Initial Register Values Consider the DDA interpolation of a parametric curve X(u) = ao + alu + a2u2 +... + a,um (Al) successive addition of the contents of the registers are accumulated in Ro. Now, suppose a total of K steps is necessary to complete the curve, then after the completion of K steps, the content of register Ro is given by K(K + 1) K(K + 1)(K + 2) (R) = a0o+Kcl+ -2! c3! + 3 2! ~~~3! IK(K + 1)(K + 2)(K + 3) 4! + K(K + 1)(K + 2)(K + 3)(K + 4) 5! a where ai is the initial value of Ri. Rearranging, we have (Ro) = ao+ a + 2+ 3+ + 5 +... }K 2 3 4 5 + 2 (1 + 2) [(1*2) + (2*3)+ (3*1)1 + {-! + a3- 4! _! _! 4! 3! 4! 51 (1 + 2 + 3) [(1*2) + (1'3) + (1*4) + (2*3) + (2*4) + (3*4)] +...}i3 4! 4 5! (1 + 2 + 3 + 4)5+...K4 5!...}K5 +... + I{... }Km (A2) By comparing (Al) and (A2), and bearing in mind the relationship that 2n in the register Ro corresponds to one unit distance (D) due to a single pulse from the DDA in the linear distance, we obtain 2n(D al (a02+a3+ a1 +a,+.... )K2 *D = a2 2!2n* 2(a3+32as+7+) * D = a3 3!2n (cr+2c+.)K4 * D = a4 (C5+,,K5 * D = as 16

Solving for al, a2, a3, and so on, we find the values to be 5!2n 5 = K D(a5 + ) 4!2n a4 = D(4 + D ~) K4 ~ D 3!2n a3 = K3D(a3+'-) 2!2n 3 a2 = K2D(a2- a3 +.) 2" a2 a3 1l = K D(alK + +*.) K*D K K2 The value of ao can be obtained by observing that the initial position is ao ideally. However, due to the fact that the position can only be in multiples of D, the actual initial position is only TRUNC(ao/D)*D, where TRUNC is a truncating operation where only the integral part of a real number is kept and the fractional part discarded. The balance of a0 which is not translated into initial-position command is stored in register Ro as ao. The value of Ca can be computed by using the following equation: ao - D*TRUNC() 2n D In the case of a rational curve, the determination of the initial values for the registers in numerator and denominator are performed separately. 17