Sage Journals: Discover world-class research

Abstract

With the outstanding geometric properties of B-spline, its recursive algorithm has been widely applied in numerical control systems to compute interpolation points and their derivatives. However, some of the B-spline basis functions in the recursive algorithm are computed repeatedly during the numerical control interpolation, which makes the interpolation efficiency greatly deteriorated. To deal with this issue, an effective approach by converting B-spline to piecewise power basis functions is proposed in this article. With the proposed power basis functions, formulas for computing interpolation point and its derivative can be established explicitly. Thus, the repeated computation due to the implicit expression can be avoided. To further improve the real-time interpolation performance, curve subdivision based on the proposed power basis functions is implemented simultaneously. Theoretical analysis shows that the computational complexity of the proposed method can be reduced by 71.1% compared with the DeBoor-Cox algorithm during the computation of interpolation points and their derivatives. Furthermore, the computational complexity of the curve subdivision can be reduced by 64.5%. Experimental results on a two-dimensional butterfly and a three-dimensional pigeon show that the proposed methods can greatly improve the computational efficiency.

Keywords

B-spline curve real-time interpolation piecewise power basis functions curve conversion curve subdivision

Introduction

Due to advantages of the simple structure and calculation, B-spline has been widely applied to describe complex curves and tool paths in the field of numerical control (NC). To realize the machining of one B-spline curve in real time, a series of operations, such as G-code interpretation, trajectory planning, acceleration/deceleration (ACC/DEC) scheduling,¹ and parametric interpolation,² are needed to be completed in a limited time period. To implement the operations mentioned above, NC systems typically adopt a multi-level structure,³ which includes NC interpreter level, NC process level, and motion control level. For example, NC process level accomplishes the complicated computation of the B-spline curve. Sometimes, the computation load may be too high to be completed in the limited time period, which affects the real-time performance of NC systems. Thus, many researchers have focused on improving the computational efficiency in NC process level including B-spline computation and interpolation process.

Considering B-spline computation complexity, for instance, the DeBoor-Cox recursive algorithm is utilized to determine interpolation points of a B-spline curve.⁴ However, when this algorithm is adopted to compute the derivatives directly,⁵ it is inevitable to compute the B-spline basis functions repeatedly, which leads to redundant computation. Grabowski and Li⁶ proposed an iterative formula to implement the DeBoor-Cox formula. Qin⁷ suggested an explicit matrix formula and pointed out that the computational efficiency of matrix formula is higher than that of the iterative formula. Biagiotti and Melchiorri⁸ put forward an algorithm, which converts the B-spline to the piecewise polynomial. Nevertheless, the polynomial coefficients need to be determined in advance according to the derivatives of B-spline nodes. Sun et al.⁹ suggested a vector extending operation to avoid the calculation of the polynomial coefficients, in which the interpolation point and the derivatives should be evaluated simultaneously. Recently, Karimi and Nategh¹⁰ proposed an algorithm to convert the B-spline into Bézier curves, but the algorithm requires the B-spline with equally spaced knots, which limits its applications.

Some researchers focus on the simplification of the interpolation process, including entire curve scanning to obtain sharp corners,¹¹ curve subdividing into several ACC/DEC schedule units (ASUs), arc length calculating, and the ACC/DEC scheduling of each ASU.¹² For example, to divide a B-spline curve quickly, Tsai et al.¹³ suggested a pre-interpolation method, which can generate a series of pre-interpolation points within the chord error. The curve is subdivided according to the curvature at these points. However, the number of the generated points varies greatly for various curves, which causes an uncertain execution time. Li et al.¹⁴ proposed a method based on the minimum splitting length. The number of points can be determined and the execution time is affordable, since the pre-interpolation points are generated with equal length interval. However, the corresponding control points need to be computed at each pre-interpolation point, which increases the complexity of computation. To obtain the velocity profile quickly, Lin et al.¹⁵ suggested a look-ahead ACC/DEC schedule. The velocity profile of ASUs could be generated step by step, but the arc length of ASU is approximated by accumulating its chord lengths, which increases the arc-length error and interpolation error. The arc-length error was minimized by the adaptive Simpson’s method,¹⁶ and the interpolation error was constrained by deducing the relationship between the interpolation point and the arc length.¹⁷ Since the curve subdivision and the arc-length calculation are two separate operations, NC system needs to implement them, respectively, which needs additional execution time.

According to the analysis above, the efficiency ofB-spline interpolation can be enhanced by two aspects: the B-spline computation and the interpolation process. To improve the algorithm efficiency, an effective approach by converting B-spline to piecewise power basis functions is proposed in this article, which is named as the curve conversion. With the proposed power basis functions, the formulas for computing interpolation point and its derivative can be established explicitly. Thus, the repeated computation due to implicit description can be avoided. To improve the interpolation process, curve subdivision based on the proposed power basis functions is implemented simultaneously. With the proposed cooperation strategy, the efficiency of real-time B-spline interpolation is greatly improved while keeping the interpolation accuracy unchanged. Theoretical analysis and experimental results show that the proposed methods can greatly improve computational efficiency.

The remainder of this article is organized as follows. The conversion of B-spline curve to piecewise power basis functions is introduced in section “Conversion from B-spline to piecewise power basis functions.” Section “Curve subdivision” shows the curve subdivision operation. Theoretical computational efficiency of both the interpolation algorithm and the interpolation process is analyzed in section “Theoretical computational efficiency comparison.” Experimental results on a 2D butterfly and a 3D pigeon are shown in section “Experiment validation,” and conclusions are given in section “Conclusion.”

Conversion from B-spline to piecewise power basis functions

A pth-order B-spline curve is defined as

\begin{matrix} C (u) = \sum_{i = 0}^{n} N_{i, p} (u) P_{i}, u \in [0, 1] \\ U = {\underset{p + 1}{\underset{︸}{0, \dots, 0}}, u_{p + 1}, \dots, u_{m - p - 1}, \underset{p + 1}{\underset{︸}{1, \dots, 1}}} \\ n = m - p - 1 \end{matrix}

(1)

where ${P_{i}}$ and $U$ denote the control points and the knot vectors, respectively, whose numbers are n + 1 and m + 1. ${N_{i, p}}$ are the pth-order B-spline basis functions defined as

\begin{matrix} N_{i, 0} = {\begin{matrix} 1 if u_{i} \leq u < u_{i + 1} \\ 0 otherwise \end{matrix} \\ N_{i, p} = \frac{u - u_{i}}{u_{i + p} - u_{i}} N_{i, p - 1} (u) + \frac{u_{i + p + 1} - u}{u_{i + p + 1} - u_{i + 1}} N_{i + 1, p - 1} (u) \end{matrix}

(2)

Substituting equation (2) into equation (1), it gets the DeBoor-Cox recursive algorithm⁴

\begin{matrix} C (u) = P_{j}^{(p)} \\ P_{i}^{(h)} = {\begin{matrix} (1 - α_{i}^{h}) P_{i - 1}^{(h - 1)} + α_{i}^{h} P_{i}^{(h - 1)} & h = 1, 2, \dots, p \\ P_{i} & h = 0 \end{matrix} \\ α_{i}^{h} = \frac{u - u_{i}}{u_{i + p + 1 - h} - u_{i}}, u \in [u_{j}, u_{j + 1}], i = j - p + h, \dots, h \end{matrix}

(3)

where $j$ is the index of $u$ in $U$ . When ${P_{i}}$ and $U$ are given, for any $u_{i} \in [0, 1]$ , interpolation point $C (u_{i})$ can be calculated by equation (3). Similarly, by utilizing the DeBoor-Cox algorithm, $C' (u_{i})$ and $C ″ (u_{i})$ , the first and the second derivatives of $C (u_{i})$ , can also be determined. Table 1 shows the statistics of arithmetic operations for DeBoor-Cox algorithm.¹⁸

Table 1.

Summary of the arithmetic operations for DeBoor-Cox algorithm.

Arithmetic operation	Interpolation point		First derivative		Second derivative
	$α_{i}^{h}$	$C (u)$	$α_{i}^{h}$	$C' (u)$	$α_{i}^{h}$	$C ″ (u)$
Addition	0	$pAD / 2$	0	$pBD / 2$	0	$BCD / 2$
Subtraction	$pA$	$pAD / 2$	$pB$	$pBD / 2 + N_{α}$	$BC$	$BCD / 2 + 2 N_{α}$
Multiplication	0	$pAD$	0	$pBD + N_{α}$	0	$BCD + 3 N_{α}$
Division	$pA / 2$	$N_{α}$	$pB / 2$	$N_{α}$	$BC / 2$	$N_{α}$

N_α: the dimension of N-space; p: the degree of B-spline curve; p = 3; A: p + 1; B: p − 1; C: p − 2; $D : N_{α} + 1$ .

With the DeBoor-Cox algorithm, NC system can implement interpolation process, including curve subdivision, arc-length calculation, and interpolation point generation. However, these operations require to call the DeBoor-Cox algorithm repeatedly, which causes redundant computation. It is noted that the B-spline curve can be converted into piecewise power basis functions as

{\begin{matrix} H_{0} (u) = a_{0, 0} u^{p} + a_{0, 1} u^{p - 1} + \cdot \cdot \cdot + a_{0, p - 1} u + a_{0, p}, u \in [0, u_{p + 1}] \\ H_{1} (u) = a_{1, 0} u^{p} + a_{1, 1} u^{p - 1} + \cdot \cdot \cdot + a_{1, p - 1} u + a_{1, p}, u \in [u_{p + 1}, u_{p + 2}] \\ \cdot \cdot \cdot \\ H_{j} (u) = a_{j, 0} u^{p} + a_{j, 1} u^{p - 1} + \cdot \cdot \cdot + a_{j, p - 1} u + a_{j, p}, u \in [u_{p + j}, u_{p + j + 1}] \\ \cdot \cdot \cdot \\ H_{m - 2 p - 1} (u) = a_{m - 2 p - 1, 0} u^{p} + a_{m - 2 p - 1, 1} u^{p - 1} + \cdot \cdot \cdot + a_{m - 2 p - 1, p - 1} u + a_{m - 2 p - 1, p}, u \in [u_{m - p - 1}, 1] \end{matrix}

(4)

Then, it becomes much easier to obtain $C (u_{i})$ . NC system usually adopts the dichotomy method to determine the interval of $u_{i}$ , and $H (u_{i})$ is achieved by substituting $u_{i}$ into its corresponding power basis function $H_{j} (u)$ . Table 2 gives the operation statistics for the power basis function algorithm.

Table 2.

Summary of the arithmetic operations for the power basis function algorithm.

Arithmetic operation	Interpolation point	First derivative	Second derivative
	$C (u)$	$C' (u)$	$C ″ (u)$
Addition	$p \times N_{α}$	$(p - 1) \times N_{α}$	$(p - 2) \times N_{α}$
Multiplication	$p (p + 1) N_{α} / 2$	$(p - 1) (p + 2) N_{α} / 2$	$(p - 2) (p + 3) N_{α} / 2 + N_{α}$

N_α: the dimension of N-space; p: the degree of B-spline curve; $p \geq 3$ .

From Tables 1 and 2, the piecewise power basis functions provide less complexity in computation and are easier to be implemented in NC system.

To clearly show the conversion procedure, an example with the cubic B-spline is described in detail below.

Step 1. Convert the B-spline curve $C (u)$ into $m - 2 p$ Bézier curves. For instance, the kth Bézier curve is expressed as

C_{k}^{Bezier} (u) = \sum_{i = 0}^{3} B_{i, 3} (u) Q_{k, i}, u \in [u_{p + k}, u_{p + k + 1}]

(5)

where ${B_{i, 3} (u)}$ are the cubic Bernstein polynomials. ${Q_{k, i}}$ are the control points of the Bézier curve, which can be easily determined by the knot intersection method.¹⁴

Step 2. Convert each Bézier curve to power basis function by matrix multiplications, for example, the kth power basis function $H_{k} (u)$ can be seen in equation (6)

\begin{matrix} H_{k} (u) = a_{k, 0} u^{3} + a_{k, 1} u^{2} + a_{k, 2} u + a_{k, 3}, u \in [u_{p + k}, u_{p + k + 1}] \\ c = \frac{1}{u_{p + k + 1} - u_{p + k}}, d = \frac{- u_{p + k}}{u_{p + k + 1} - u_{p + k}} \\ [\begin{matrix} α_{k, 0} \\ α_{k, 1} \\ α_{k, 2} \\ α_{k, 3} \end{matrix}] = [\begin{matrix} - c^{3} & 3 c^{3} & - 3 c^{3} & c^{3} \\ 3 (c^{2} - d c^{2}) & (9 d c^{2} - 6 c^{2}) & (3 c^{2} - 9 d c^{2}) & 3 d c^{2} \\ 3 (2 cd - c - c d^{2}) & (9 c d^{2} - 12 cd + 3 c) & 6 cd - 9 c d^{2} & 3 c d^{2} \\ 1 - 3 d + 3 d^{2} - d^{3} & 3 (d - 2 d^{2} + d^{3}) & 3 (d^{2} - d^{3}) & d^{3} \end{matrix}] [\begin{matrix} Q_{k, 0} \\ Q_{k, 1} \\ Q_{k, 2} \\ Q_{k, 3} \end{matrix}] \end{matrix}

(6)

Step 3. Calculate the first and second derivatives of $H_{k} (u)$ due to equation (6)

\begin{matrix} {H'}_{k} (u) = 3 α_{k, 0} u^{2} + 2 α_{k, 1} u + α_{k, 2}, u \in [u_{p + k}, u_{p + k + 1}] \\ {H ″}_{k} (u) = 6 α_{k, 0} u + 2 α_{k, 1} \end{matrix}

(7)

Through above-mentioned steps, a set of power basis functions, whose number is $m - 2 p$ , are obtained, and the statistics of arithmetic operations for the curve conversion are shown in Table 3. There is no conversion error during the conversion procedure theoretically. However, computational error does exist in a discrete system, which deteriorates the computational accuracy of the interpolation points. To deal with this problem, the error can be constrained by setting an error threshold $ε$ .

Table 3.

Summary of the arithmetic operations for curve conversion.

Process	Arithmetic operation
	Addition	Subtraction	Multiplication	Division
Step 1	$3 N_{α} (m - 2 p - 1)$	$6 (m - 2 p - 1)$	$6 N_{α} (m - 2 p - 1)$	$2 (m - 2 p - 1)$
Step 2	$(4 p N_{α} - N_{α} + 3) \times (m - 2 p)$	$(4 p N_{α} - N_{α} + 3) \times (m - 2 p)$	$(4 p N_{α} + 4 N_{α} + 29) \times (m - 2 p)$	$2 \times (m - 2 p)$

N_α: the dimension of N-space; p: the degree of B-spline curve; p = 3.

Curve subdivision

As shown in Figure 1, the traditional interpolation process first executes the pre-interpolation for scanning curvature extrema. Then, the curve subdivision utilizes these extrema to divide the curve and generate ASUs. After that, the arc-length calculation and the ACC/DEC schedule are implemented for each ASU to generate the corresponding velocity profile. In the process, computational load depends on operating the pre-interpolation and calculating arc lengths.¹⁷ With the power basis functions, we can integrate curve subdivision and arc-length calculation to improve the computational efficiency. The process is described in detail as below:

Step 1. Adopt the adaptive Gauss-Legendre quadrature rule to calculate the arc length of $C (u)$

Δ s = \int_{- 1}^{1} f (x) dx \approx \sum_{i = 0}^{5} A_{i} f (x_{i})

(8)

Equation (8) denotes the six-point Gauss-Legendre quadrature rule, and $f (x)$ is an integrand. $x$ is an integral variable over an interval $x \in [- 1, 1]$ , and ${x_{i}}$ are Gauss nodes. ${A_{i}}$ are weights, and $Δ s$ is the analytic result. The calculation of the arc length of $C (u)$ is as follows: when iterative process starts, the interval is subdivided equally and equation (8) is applied to each subintervals. The approximation value is then obtained by summing the recursive subintervals. In this way, the iterative process is continued until $δ$ (the absolute value of the difference between iterations) is less than the required tolerance $τ$ . After the iterative process, the arc length of $C (u_{i})$ is calculated. Meanwhile, $n$ subintervals and $6 n$ relevant parameters are also obtained. Figure 2(a) shows the schematic diagram of subintervals, and one of the subintervals, named as $SU B_{k}$ , is illustrated in Figure 2(b). According to equation (8), the arc length in the $SU B_{k}$ is expressed as

\begin{matrix} L_{k} = \int_{a}^{b} ‖ C' (u) ‖ du \approx 0.5 (b - a) \sum_{i = 0}^{5} A_{i} ‖ H' (u_{i}) ‖ \\ u_{i} = 0.5 (b - a) x_{i} + 0.5 (b + a), i = 0, \dots, 5 \end{matrix}

(9)

where $a$ and $b$ are the lower and the upper limit of integral. $b - a$ is the integral length, which is equal to 1/n. ${u_{i}}$ are the mapping values of Gauss nodes ${x_{i}}$ .

Figure 1.

Comparison of the traditional interpolation process and the proposed interpolation process.

Figure 2.

The schematic diagram of SUBs: (a) the distribution of SUBs and (b) the enlarged view of SUBk.

Step 2. Determine curvature extrema and inflection points from the set of generated parameters.

$H' (u_{i})$ has been acquired by equation (9). NC system only needs to calculate $H ″ (u_{i})$ by equation (7), the corresponding curvature $ρ (u_{i})$ will be rapidly obtained

ρ (u_{i}) = \frac{‖ H' (u_{i}) \times H ″ (u_{i}) ‖}{{‖ H' (u_{i}) ‖}^{3}}

(10)

Through scanning the curvatures of these 6n parameters, the curvature extrema of the curve can be roughly determined. Obviously, the tolerance $τ$ is smaller, the number of SUBs gets greater, and the accuracy of curvature extreme also achieves higher. However, the computational load is greatly increased. To balance the accuracy and the performance, some parameters are selected from these 6n parameters in this article. As long as the distribution of interpolation points of the selected parameters is as even as possible, an effective curvature scanning with less calculation will be realized. The distribution of interpolation points can be roughly evaluated by the Gauss node increment, and the proof is given below.

Define the parameter increment, the arc length increment, and the Gauss node increment as $Δ u = u_{i} - u_{i - 1}$ , $Δ L = ‖ C (u_{i + 1}) - C (u_{i}) ‖$ , and $Δ x = x_{i} - x_{i - 1}$ , respectively. Assume $u_{i}, u_{i - 1} \in [u_{p + k}, u_{p + k + 1}]$ , when $Δ u << 1$ , it obtains $Δ L \approx ‖ a_{k, 2} ‖ Δ u$ by equation (6). In addition, according to equation (9), it derives $Δ u = Δ x / 2 n$ . Thus, when $Δ u << 1$ , the arc length increment is $Δ L \approx (‖ a_{k, 2} ‖ / 2 n) Δ x = q Δ x$ . That is, in each subinterval, Gauss node increment can approximately reflect the distribution of $C (u_{i})$ . After that, based on the Gauss node increment, the principle of the parameter selection is as follows:

The parameters at ends of the curve (i.e. 0 and 1).

The parameters $u_{1}$ , $u_{2}$ , $u_{3}$ , and $u_{4}$ in each SUB.

According to this principle, the distribution of interpolation points of selected parameters is much more even, which can be seen in Figure 3(a). Meanwhile, the distribution of that is approximated as follows: 0.1694L₀, 0.2113L₀, 0.2386L₀, 0.2113L₀, 0.1694(L₀ + L₁), 0.2113L₁, 0.2386L₁, 0.2113L₁, …, 0.1694(L_n_− 2 + L_n_− 1), 0.2113L_n_− 1, 0.2386L_n_− 1, 0.2113L_n_− 1, and 0.1694L_n_− 1.

Figure 3.

The distribution of the selected points: (a) the selected points on the curve and (b) the selected points on the SUB_k.

In this way, NC system only calculates 4n + 2 curvatures in total. Compared with the traditional operation that only curvature extrema are scanned, the proposed operation will scan both curvature extrema and inflection points (as seen in Figure 4). Thus, when ASUs are generated based on these points, only four forms of velocity profiles are needed (i.e. acceleration phase; acceleration and constant velocity phase; deceleration phase; and deceleration and constant phase), which greatly simplifies the ACC/DEC schedule.

Figure 4.

Different operations for curve subdivision: (a) the traditional operation and (b) the proposed operation.

Step 3. Calculate the arc length of ASU.

First, determine the subinterval where the endpoint of ASU is located. Then, calculate the relative arc length of ASU in this subinterval by equation (9). After that, accumulate all arc lengths of the previous subintervals. Therefore, the arc length of ASU is estimated with lower computational complexity.

Theoretical computational efficiency comparison

In this section, a cubic B-spline curve is used in the comparison, where the degree of the curve p = 3, the number of control points ${P_{i}}$ (i = 0, …, n) and knots ${u_{i}}$ (i = 0, …, m) are n = 9 and m = 13, respectively. Meanwhile, to quantify the comparison of the computational efficiency, the instruction cycle (latency) is used as the evaluation index. Note that the latencies for addition/subtraction, multiplication, division, and square root are 3, 5, 24, and 34, respectively:¹⁹

Computational efficiency of curve conversion: in section “Conversion from B-spline to piecewise power basis functions,” the arithmetic operations for the curve conversion have been summarized (as seen in Table 3). Thus, the latencies of curve conversion operation can be easily determined, which is equal to 5242. After the curve conversion operation is completed, the B-spline curve is converted to the piecewise power basis functions. With the power basis functions, $C (u)$ , $C' (u)$ , and $C ″ (u)$ can be quickly computed, and their latencies are 117, 93, and 69, respectively. However, with the DeBoor-Cox algorithm,⁴ the latencies of calculating $C (u)$ , $C' (u)$ , and $C ″ (u)$ are 456, 288, and 233, respectively. Compared with the DeBoor-Cox algorithm, the latencies of the proposed curve conversion algorithm can be reduced by about 71.1% on average (as seen in Table 4). Although there is a 5242 latency for curve conversion before the computation of interpolation points and their derivatives, we may need to compute the points and their derivative for up to several hundred or thousand times. Moreover, only one curve conversion is needed for a curve. Thus, compared to the total time, the curve conversion is acceptable for real applications.

Computational efficiency of curve subdivision: since the arc length is calculated by the adaptive Gauss–Legendre quadrature rule during curve subdivision operation, some parameters of the quadrature formula should be configured before iteration. Here, the initial number of intervals and the maximal number of intervals are 7 and 10, respectively. Meanwhile, assume that the curve is divided into three ASUs after iteration. Accordingly, the latency of the curve subdivison operation based on the power basis functions is 46,515. Under the same parameters, the latency of pre-interpolation algorithm¹⁷ is 130,953. Thus, the latencies of the proposed curve subdivision operation can be reduced by 64.5% compared with that of pre-interpolation algorithm.

Table 4.

Comparison of the latencies between the proposed algorithm and the DeBoor-Cox algorithm.

Algorithm	$C (u)$	$C' (u)$	$C ″ (u)$
Proposed algorithm (latency)	117	93	69
DeBoor-Cox algorithm⁴ (latency)	456	288	233

Experiment validation

Platform construction

In this section, a motion control platform is set up to verify the proposed algorithm, which consists of an upper computer, a motion controller, and a motion platform:

An upper computer: the computer is configured with Core i5-3470 3.2 GHz CPU, Kingston^® DDR3 1600 MHz 8GB memory, and a PC-based NC system. The main parts of NC system, including a smoothing module, a look-ahead ACC/DEC schedule module, an interpolation module, and a timer module, are coded with Visual Studio 2010. Here, the cycle of NC process level is set to 10 ms based on the Win 32 multimedia timer²⁰ (with 1-ms resolution) and the theoretical interpolation cycle is 1 ms. Thus, NC system can generate 10 interpolation points in each interrupt service routine (ISR).

To generate the interpolation points in real time, all modules are required to collaborate with each other. This collaboration mechanism is defined as the real-time interpolation strategy. Meanwhile, a sliding window is introduced here, which can cover the maximum number of $k$ line segments. With the aid of sliding window, the real-time interpolation strategy is implemented as follows:

(a) In the first cycle (expressed as $I_{0}$ in Figure 5(b)), NC system puts $k$ line segments into the window (e.g. window 0 in Figure 5(a)). Then, the timer module activates ISR and enables other modules in ISR. First, the smoothing module smoothes these line segments and generates a cubic B-spline curve $C (u)$ (the solid red curve in window 0), which will be converted into power basis functions. Second, the curve subdivision operation is performed, which generates ASUs and executes the arc-length calculation. Then, the look-ahead ACC/DEC module obtains the velocity profiles. At last, based on ASUs and feedrate profiles, the interpolation module generates the specified number of interpolation points (10 points) and utilizes the points to guide the tool movement along the curve $C (u)$ .

(b) In the next cycle, NC system examines whether the tool has reached the window-sliding position (the position where the window 0 and windows 1 overlap). If not reached, only the interpolation module works in this ISR (the $I_{1}$ in Figure 5(b)). Accordingly, new interpolation points, with the same number as before (10 points), are generated and make the tool continue to move. Otherwise, the window slides forward (the window 1 in Figure 5(a)). That is, the line segments in front of the window-sliding position are removed from the window, and the same number of new line segments are added into the window. Then, all modules are involved.

A motion controller: the controller is used to implement the fine interpolation and the closed-loop control, which is based on TMS 320F2812 DSP and Cyclone II FPGA.

A motion platform: the platform consists of an AC servo driver unit and a X-Y platform. The AC servo driver unit consists of two YASKAWA SGDV-2R8A01A Series AC servo drivers and two SGMJV-04ADA21 Series AC servo motors. The movement range of the X-Y platform is 200 mm × 200 mm.

Figure 5.

The schematic diagram of real-time interpolation: (a) the machining curve and (b) the program of the NC system.

Interpolation experiment and efficiency analysis

Figure 6 shows a 2D butterfly and a 3D pigeon, which compose of 501 and 270 line segments, respectively. The real-time interpolation for them will be implemented here to verify the effectiveness of the proposed method. In the experiments, the NC system with the proposed method is named as “NC system 1” and the NC system without the proposed curve conversion method is called as “NC system 2.” As seen in Figure 7, after the interpolations for the butterfly pattern, NC system 1 and NC system 2 generate the same number of interpolation points (the number is 34,082). The actual velocity curves are illustrated in Figure 8(a). As shown in the figure, they are completely overlapped. Similarly, for the pigeon pattern, two NC systems have the same generated points (the number is 43,632), the actual velocity curves of which are also completely overlapped as shown in Figure 8(b).

Figure 6.

Instances of tool path curves: (a) 2D butterfly and (b) 3D pigeon.

Figure 7.

The experiments of real-time interpolation for 2D butterfly pattern.

Figure 8.

The velocity curves after real-time interpolation: (a) the butterfly pattern and (b) the pigeon pattern.

To exactly record the code execution time in the experiments, two Win 32 API functions, QueryPer $formanceFrequency ()'$ and $QueryPerformance Counter ()'$ are used.²⁰

For the butterfly pattern, the execution time of the curve conversion and curve subdivision with NC system 1 is illustrated in Figure 9(a). The execution time of the curve subdivision with NC system 2, which directly utilizes the DeBoor-Cox algorithm,⁵ is exhibited in Figure 9(b). Clearly, in Figure 9, (a) and (b) have the same shape but different execution times, especially at the moments when it is required to frequently calculate interpolation points and derivatives. As seen in Figure 9(a), arrow 1 denotes the execution time of curve conversion and curve subdivision in the first cycle. According to the sliding window characteristics, $k$ line segments are put into the empty window in this cycle. Then, all these line segments should be smoothed. In other cycles, only new line segments need to be smoothed. Thus, the execution time of curve generation as well as curve subdivision in the first cycle is larger than that in any other cycle. Arrow 2 represents the execution times of curve conversion and curve subdivision in three cycles, which need 17, 22, and 25 iterations to complete the arc-length calculation, respectively. Nevertheless, in other cycles, the number of iterations is less than 10 times. This phenomenon also appears in NC system 2 (arrows in Figure 9(b)). Therefore, although the B-spline curve and piecewise power basis functions can be transformed from one to another, their execution efficiencies are significantly different.

Figure 9.

Execution times of curve conversion and curve subdivision with NC system 1 and execution times of curve subdivision with NC system 2: (a) interpolation for butterfly pattern by NC system 1, (b) interpolation for butterfly pattern by NC system 2,(c) interpolation for pigeon pattern by NC system 1, and (d) interpolation for pigeon pattern by NC system 2.

The traditional interpolation process needs to implement the pre-interpolation, and the computation is obviously much more than the proposed process. In this article, we simplify the traditional interpolation process by combing the curve subdivision with arc-length calculations, so the execution time with NC system 1 is less than 10 ms, as shown in Figure 9(a). Here, both the error threshold $ε$ and the tolerance $τ$ are set to $10^{- 5} mm$ . Figure 10(a) shows the distribution of the generated ASUs. In the figure, each ASU is a curve segment from the first symbol “+” to the next one.

Figure 10.

The distribution of the ASUs: (a) on the butterfly pattern and (b) on the pigeon pattern.

A comprehensive evaluation for the proposed method is also given by recording the execution time of each ISR. When real-time interpolations for butterfly pattern is finished, NC system 1 and NC system 2 all generate 3408 ISRs. According to the previous analysis, some execution times of ISRs may exceed the time interval (i.e. 10 ms). In this case, the ISR will not quit until all calculations are completed. Obviously, it will impact the machining quality and efficiency. As shown in Figure 11(a), when NC system 1 is applied, all ISRs are less than the time interval. However, NC system 2 exists the problem that the execution time of some ISRs exceeds the time interval which is 10 ms, as shown in Figure 11(b). From the previous analysis, the execution time with NC system 2 is 10 times higher than that with NC system 1.

Figure 11.

Different execution times of the real-time interpolation: (a) interpolation for butterfly pattern by NC system 1, (b) interpolation for butterfly pattern by NC system 2, (c) interpolation for pigeon pattern by NC system 1, and (d) interpolationfor pigeon pattern by NC system 2.

Similarly, for the pigeon pattern, Figure 9(c) and (d) denotes the execution times of the curve conversion and curve subdivision with NC system 1 and with NC system 2, respectively. Figure 10(b) gives the distribution of the generated ASUs. In addition, the execution times of each ISR with NC system 1 and NC system 2 are illustrated in Figure 11(c) and (d), respectively. From these figures, it can be seen that the efficiency of the NC system can be greatly improved by the proposed method.

Conclusion

This article proposes an algorithm for real-time interpolation based on piecewise power basis functions, which converts the B-spline into power basis functions. Since the interpolation point and the derivative formula can be established explicitly based on the power basis functions, interpolation points and their derivatives can be rapidly computed. Meanwhile, the curve subdivision based on the piecewise power basis functions is implemented, which can further simplify the interpolation process. Thus, the real-time performance for the interpolation can be improved greatly. Theoretical analysis illustrates that the latencies of the proposed curve conversion method and the curve subdivision method can be reduced by 71.1% and by 64.5% compared with the DeBoor-Cox algorithm and the pre-interpolation method. Experiments are carried out with a 2D butterfly and a 3D pigeon, and the results show that the proposed method can greatly improve the computational efficiency.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported, in part, by National Key Basic Research Program of China (Grant No. 2013CB035804), National Science and Technology Major Project (Grant No. 2015ZX04005006), and State Key Lab of Digital Manufacturing Equipment & Technology (Grant No. DMETKF2016002).

ORCID iD

Hui Wang

References

Huang

. Parallel acceleration/deceleration feedrate scheduling for computer numerical control machine tools based on bi-directional scanning technique. Proc IMechE, Part B: J Engineering Manufacture. Epub ahead of print 24 May 2017. DOI: 10.1177/0954405417706997.

Zhang

. Fast parametric curve interpolation with minimal feedrate fluctuation by cubic B-spline. Proc IMechE, Part B: J Engineering Manufacture 2016; 24(1): 1–11.

Wang

Yau

. Universal real-time NURBS interpolator on a PC-based controller. Int J Adv Manuf Tech 2014; 71(1–4): 497–507.

Jaen-Cuellar

De Santiago-Perez

Romero-Troncoso

. A hardware-software system for coordinated multi-axis control based on a non uniform rational B-splines interpolator applied to industrial computer numerically controlled machines. Proc IMechE, Part I: J Systems & Control Engineering 2012; 226(6): 831–840.

Chen

Liao

Cheng

. Design and implementation of integrated non-uniform rational B-spline and digital differential analyser interpolators for computerized numerical control servocontrollers. Proc IMechE, Part C: J Mechanical Engineering Science 2007; 221(9): 1075–1087.

Grabowski

. Coefficient formula and matrix of nonuniform B-spline functions. Comput Aided Design 1992; 24(12): 637–642.

Qin

. General matrix representations for B-splines. Visual Comput 2000; 16(3–4): 177–186.

Biagiotti

Melchiorri

. Trajectory planning for automatic machines and robots. Berlin: Springer, 2008.

Sun

Chai

Fan

. A fast algorithm to compute B-splines and its application in curve fitting. Key Eng Mat 2008; 392–394: 743–749.

10.

Karimi

Nategh

. Development of a novel adaptive nonuniform rational basis spline interpolator with limited kinematic error for hexapod machine tools. Proc IMechE, Part B: J Engineering Manufacture 2013; 228(228): 319–327.

11.

Yong

Narayanaswami

. A parametric interpolator with confined chord errors, acceleration and deceleration for NC machining. Comput Aided Design 2003; 35(13): 1249–1259.

12.

Lee

Lin

Pan

. The feedrate scheduling of NURBS interpolator for CNC machine tools. Comput Aided Design 2011; 43(6): 612–628.

13.

Tsai

Nien

Yau

. Development of an integrated look-ahead dynamics-based NURBS interpolator for high precision machinery. Comput Aided Design 2008; 40(5): 554–566.

14.

Qiu

Zhang

. A practical real-time non-uniform rational B-spline curve interpolator for computer numerical control machining. Proc IMechE, Part C: J Mechanical Engineering Science 2012; 226(4): 1068–1083.

15.

Lin

Tsai

Yau

. Development of real-time look-ahead algorithm for NURBS interpolator with consideration of servo dynamics. In: Proceedings of the 46th IEEE conference on decision and control, vol. 12, New Orleans, LA, 12–14 December 2007, pp.1862–1867. New York: IEEE.

16.

Heng

Erkorkmaz

. Design of a NURBS interpolator with minimal feed fluctuation and continuous feed modulation capability. Int J Mach Tool Manu 2010; 50(3): 281–293.

17.

Liu

Huang

Yin

. Development and implementation of a NURBS interpolator with smooth feedrate scheduling for CNC machine tools. Int J Mach Tool Manu 2014; 87: 1–15.

18.

Zhao

Zhu

Xiong

. Development of FPGA based NURBS interpolator and motion controller with multiprocessor technique. Chin J Mech Eng 2013; 26(5): 940–947.

19.

Intel. Intel 64 and IA-32 architectures optimization reference manual: system programming guide. Santa Clara, CA: Intel, 2011.

20.

Kuperberg

Krogmann

Reussner

. TimerMeter: quantifying properties of software timers for system analysis. In: Proceedings of the sixth international conference on the quantitative evaluation of systems, Budapest, 13–16 September 2009, pp.85–94. New York: IEEE.

Research of the real-time interpolation based on piecewise power basis functions