Global toolpath smoothing for CNC machining based on B-spline approximation with tool tip position adjustment

B-spline curve approximation is widely used for fitting linear toolpaths to improve machining quality and efficiency in CNC machining. For a high-quality smoothing method, the control of both approximation error and curve curvature needs to be taken into account. To reduce the maximum curve curvature while meeting the precision requirements, a B-spline approximation scheme with tool tip position adjustment is proposed in this paper. Toolpaths are first divided into several subdivisions according to the discrete curvature of each tool tip point. For subdivisions that need adjustment, tool tip positions are adjusted to minimize the maximum discrete curvature. An existing approximation method named energy-term-incorporated progressive and iterative approximation for least square fitting (ELSPIA) is selected and improved to fit the adjusted toolpaths and lower the chord errors. For improving the numerical stability of the ELSPIA method, the way of searching appropriate foot point parameters is discussed in this paper. Deviation expansion factors of tooltip points are also defined to make the ELSPIA method suitable for fitting the adjusted toolpaths. Both simulations and experimental studies are conducted to prove that the proposed method can significantly decrease the maximum curvature of fitted B-spline curves and improve the machining efficiency without exceeding the tolerance. For example, experimental results show that for the tested butterfly and rabbit toolpaths, the proposed method can improve the machining efficiency by 2.61% and 2.37%, compared with the ELSPIA method without tool tip position adjustment.


Introduction
Nowadays, the linear toolpath (G01 blocks), composed of a set of short line segments, is the most common toolpath format used in CNC machining [1]. However, the G1 and G2 discontinuities at the junctions of two linear segments will inevitably induce frequent feedrate fluctuation and acceleration variation [2], which has a significant negative influence on machining quality and efficiency in high-speed machining. Therefore, toolpath smoothing operation is essential in numerical control systems to improve toolpath continuity [3].
Toolpath smoothing methods can be classified into two types: local transition and global fitting methods. The local transition strategy, also called the corner smoothing strategy, replaces the sharp corners with transition segments, e.g., parametric curves. Zhao et al. [4] performed corner rounding for 3-axis toolpaths by using the cubic B-spline with 5 control points and derived the analytical expressions of the curvature extrema and approximation error. A realtime look-ahead interpolation scheme was further developed to acquire feedrate and acceleration with smooth profiles. Bi et al. [5] proposed a dual-Bézier transition algorithm to smooth the segment joints of both the translational path and rotational path in 5-axis machining. Huang et al. [6] analytically controlled the tool orientation smooth errors and parameter synchronization between tool positions and orientations of 5-axis toolpaths by utilizing a pair of cubic B-splines as transition curves. Nevertheless, the local transition methods are not suitable for manufacturing complicated freeform surfaces [7]. According to the predefined machining tolerance, the linear segment length of a complicated surface may fall into a submillimeter or even micrometer magnitude [8]. Since the curvature extremum of a transition curve is usually inversely proportional to its length [9] and the nominal velocity in feedrate scheduling is constrained by the curvature extremum [10], micro-line segments can lead to tiny transition curves and great curvature extrema. Consequently, feedrate fluctuation will be introduced during interpolation, which will weaken the effects of local transition on machining efficiency and quality [11].
The global fitting strategy, converting a set of consecutive linear segments into a smooth parametric curve, can be further divided into two kinds: interpolation and approximation methods. The B-spline curve is the mostly used parametric curve for global fitting. The interpolation methods aim at constructing a parametric curve whose control point number equals that of tool tip points, so that the curve is guaranteed to pass through each position accurately [12,13]. Therefore, two shortcomings of the interpolation methods emerge innately. Firstly, the interpolation methods are not suitable for smoothing toolpaths consisting of numerous micro-line segments. Secondly, the shape of the curve between adjacent tool tip points is hard to control, so the curve may wiggle frequently and bring out excessive curvature extremum points. By comparison, the approximation methods enable to obtain a resultant B-spline curve with fewer control points. Thus, the approximation methods are more competitive than the local transition methods and the interpolation methods to acquire ideal smoothed toolpaths for machining complicated surfaces.
The conventional B-spline fitting method, the least square fitting method, minimizes the least square fitting error or energy functions to solve the control points of the fitted B-spline curve [14]. Therefore, it is required to solve a large system of linear equations and inverse matrices, which raises challenges of numerical stability and computational load. To overcome this drawback, Lin et al. [15] proposed the progressive and iteration approximation (PIA) method for B-spline fitting. The PIA method constructs the initial B-spline curve by selecting all the data points as initial control points, then obtaining the fitted curve within tolerance by adjusting the control points iteratively and explicitly. However, the PIA method still maintains the same number of control points as the data points. To improve the PIA method, Lin and Zhang [16] developed an extended PIA format, which needs fewer control points than data points for curve fitting. After that, Deng and Lin [17] developed a new progressive and iterative approximation for least square fitting (LSPIA). For the LSPIA method, the initial control points can be arbitrarily selected for given data points, and the knot vector and control points can be adjusted during iterations with an incremental algorithm. In addition, the limit curve constructed through LSPIA is the result of the least square fitting to the same data points. Although the LSPIA method has presented its flexibility, efficiency, and robustness in handling large-size point sets [17], it fails to settle the common problem of global fitting methods: the chord errors between the linear segments and the fitted curve are not considered [18,19].
Recently, several methods have been reported to constrain chord errors. Liang et al. [20] inserted knots iteratively during the least square fitting procedure until the overall chord error was conformed. Bi et al. [21] proposed a B-spline fitting scheme taking advantage of the B-spline's strong convex hull property. The scheme splits a roughly fitted B-spline into curve segments based on the PIA method and knot insertion algorithm. Then the knot vectors and control points of the curve segments, whose chord errors exceed the tolerance, will be refined analytically. However, the final fitted B-spline includes far more control points than the given data points. Du et al. [22] established a geometric deviation model between the linear segments and the fitted B-spline curve. The overall chord error is constrained by optimizing the roughly fitted B-spline curve considering geometric deviation and smoothness. He et al. [23] proposed the ELSPIA method for B-spline fitting and chord error refinement. The ELSPIA method inherits the advantages of the LSPIA method while confining the chord error and enhancing the stability.
Moreover, dominant point selection plays an essential role in global approximation methods because appropriate dominant points can not only preserve the toolpath shape, but also significantly reduce the number of control points and accelerate the convergence. Park et al. [24,25] selected dominant points from the local curvature maximum and minimum points and generated an error-bounded B-spline curve with fewer control points. However, for this method, new dominant points should be added iteratively according to the B-spline shape information after each least square fitting step, which is time-consuming. Zhao et al. [26] chose dominant points from the local curvature maximum points and added new dominant points in flat areas according to the chordal deviation. Du et al. [22] proposed a dominant point selection scheme based on chord error, estimated by a third-order osculating helix [27] and a bi-chord error test. He et al. [23] proposed a fast and convenient dominant point selection approach for ELSPIA based on the local curvature difference maximum.
The approximation methods mentioned above pay much attention to fitting error control while the curvatures are discussed little. Du and Wang [28] reduce the curvature of transition NURBS curves by determining the weights of control points with an osculating circle. Their study demonstrates that the curvature and smoothness of transition curves can influence the machining efficiency and should not be underestimated. Our previous research [29] made full use of the tolerance band to reduce the maximum curvature of a B-spline curve with 7 control points by optimizing the positions of two control points, but it was used for local transition. To reduce the maximum curve curvature while meeting the precision requirements, a global toolpath smoothing scheme based on B-spline approximation with tool tip position adjustment is proposed in this paper. As shown in Fig. 1, specific tool tip points will be identified and adjusted from 3-axis toolpaths based on the discrete curvature information. Then the ELSPIA method will be used to fit the adjusted tool tip points and control the chord errors. The output B-spline curve can be applied in CNC interpolation and machining, which will not be discussed in this paper.
The rest part of this paper is organized as follows: in Section 2, the algorithm for tool tip position adjustment is proposed. The B-spline fitting scheme, including dominant points selection, improved ELSPIA method, and approximation error estimation are introduced in Section 3. In Section 4, experiments are conducted to validate the effectiveness of the proposed method. Conclusions and further discussions are included in Section 5.

Preprocess operation before B-spline fitting
The proposed global smoothing methodology contains two stages: the preprocess operation and the B-spline curve fitting. To decrease the maximum curvature of the fitted B-spline curve, the preprocess operation is conducted by adjusting the given tool tip positions based on the discrete curvatures of tool tip points.

Curvature estimation for discrete tool tip points
For a regular curve (s) , where s is the arc-length parameter, the curvature of (s) at s 0 is defined as Suppose the arc-length parameters corresponding to points j n j=0 are s j n j=0 . Then, the second derivative of the curve (s) at s j can be approximately obtained by computing the second difference quotients as follows: The relationship between the second derivative and the second difference quotient is So that the discrete curvature of the point j can be computed by where j is the angle between vector ����������� ⃗ j−1 j and vector ����������� ⃗ j j+1 .

Selection scheme for tool tip positions that need adjustment
Since the number of tool tip points may be enormous, it is time-consuming to adjust every tool tip position. A scheme is proposed to judge whether a tool tip position needs to be adjusted. Firstly, the tool tip points are classified into several subdivisions by split points. Then the subdivisions that need adjustment are selected by comparing their maximum curvature with a predefined threshold. The endpoints, 0 and n , are automatically considered split points, and other split points are chosen from the local curvature minimum points. The curvature of each tool tip point is evaluated by Eq. (4) to obtain j n j=0 , the sequence of discrete curvature values. Since the trajectory discretization and curvature evaluation error may result in the fluctuation of discrete curvature, the split points should be selected with noisy points filtering [26]. Therefore, the point j (0 < j < n) will be selected as a split point if the following conditions are satisfied, and the diagram of the split point selection scheme is shown in Fig. 2.
where j is the discrete curvature of the point j , l is a predefined curvature upper bound of potential split points, p j and n j denote the previous and the next local curvature maximum value of j in the sequence j n j=0 , and f is the filtering parameter.
Suppose that a subdivision, determined by split points s and e , consists of tool tip points j e j=s , 0 ≤ s < e ≤ n , it will be selected for adjustment if the curvature maximum of its tool tip points exceeds a predefined curvature threshold c .
If a subdivision is selected for further adjustment, all the points in this subdivision, except for two split points, need to be adjusted. To illustrate this problem more intuitively, a case study is conducted and shown in Fig. 3. Originally, 702 tool tip points are sampled from a butterfly curve. The filtering parameter f and the curvature upper bound l are both 0.2 mm −1 , and the curvature threshold c is 2 mm −1 . After applying the selection algorithm, the curve is divided into 28 subdivisions by 29 split points, and 112 points from 8 subdivisions need to be adjusted.

Tool tip position adjustment
Suppose the tool tip points in the subdivision j e j=s , 0 ≤ s < e ≤ n , need to be adjusted. s and e are the split points. Tool tip position adjustment aims at decreasing the curvature maximum of the subdivision. An optimization model, which is shown in Fig. 4a, is developed as follows: Subdivision needs adjustment Subdivision no need for adjustment Split points , j is the optimization variable denoting the translation vector of each adjusted point. During the optimization process, two split points remain in the original position, and d stands for the magnitude upper bound of j . Generally, d should be set smaller than the tolerance of the B-spline curve fitting. Furthermore, ′ j can be calculated by Eq. (8) based on Eq. (4).
By solving the optimization problem for translation vectors j using the sequential quadratic programming method, the new tool tip positions can be obtained. A diagram of the fitted curves of original and adjusted points, constrained by a fixed tolerance band (or an envelope volume for a 3D curve), is shown in Fig. 4b. The fitted curve of the adjusted tool tip points is expected to have a smaller maximum curvature, compared with the fitted curve of the original tool tip points.
For the butterf ly curve, the tool tip positions before and after adjustment are shown in Fig. 5a, with d = 0.02 mm. Additionally, it is revealed in Fig. 5b that the local discrete curvature maximums of adjusted tool tip points decrease significantly, with a maximal reduction of 64.2% from 31.8 to 11.4.

Preliminaries of the B-spline curve
Note that the cubic B-spline curve is used for curve fitting in this paper. A p-degree B-spline curve, with m + 1 control points 0 i m i=0 ∈ ℝ 3 and a non-periodic knot vector U = u 0 , u 1 , ⋯ , u m+p+1 , is defined as are the basic functions of the p-degree B-spline curve defined on U, and can be recursively defined as According to Ref. [14], the curvature of a B-spline curve at parameter t can be calculated by If there exist several foot points, then the nearest foot point j to j will be selected as the ideal one, and the pointto-curve distance of j will be defined as

Appropriate dominant points selection
B-spline fitting with appropriate dominant points can not only preserve the toolpath shape, but also significantly reduce the number of control points and accelerate the convergence. In this paper, the local difference curvature maximum method [23] is adopted for dominant points selection. The difference curvature of the tool tip point j is defined as  where j is the discrete curvature of j , and s is a positive constant defining a neighborhood. A tool tip point is selected as a dominant point when the following condition is satisfied.
To obtain a uniform distribution of the dominant points, new points will be added if there are more than 2 s points between two adjacent dominant points. Assume a and b are two dominant points with b − a ≥ 2s , then the newly added dominant point is c where c = round a+b 2 .

ELSPIA-based B-spline fitting method
The ELSPIA method [23] presents an intuitive curvefitting way to avoid the numerical instability problem of conventional least square fitting methods when solving matrix inversion. Furthermore, it can not only adjust the point parameters, knot vector, and control points during (18) j > 0 and j > j−1 and j > j+1 the fitting process, but also control the chord errors. Therefore, the ELSPIA method is selected and utilized for B-spline curve fitting in this paper.
Before curve fitting, the knot vector and parameters of all the adjusted tool tip points should be determined. The parameters t j n j=0 of the adjusted tool tip points j n j=0 can be obtained through the normalized accumulated chord parameterization method [14].
The initial control points 0 i m i=0 are assumed to be the dominant points selected in Section 3.2, then averaging technique [25] can be used to determine the knot vector.  j,p (t)dt , which can be estimated by numerical integration. is the weight of energy term and can be evaluated as = trace( T ) trace( ) . To reduce the computational load, the procedure of ELSPIA is divided into two modules: B-spline fitting and chord error refinement. During B-spline fitting, only the least square fitting error term ( . A theoretical best for the iteration to reach its fastest convergence rate is = 2 min + max [17], where max and min are the largest and smallest eigenvalues of the matrix T respectively. Conventionally, the iteration terminates when the iteration number reaches a predefined value or the precision requirement is satisfied.
where tol is the fitting tolerance. However, the resultant curve differs slightly from the result of fitting the original points if fitting the adjusted tool tip points with a constant tolerance. This is because d, the magnitude upper bound of translation vectors mentioned in Section 2.3, is less than the tolerance tol . To significantly reduce the maximum curvature of the fitted curve, the fitting tolerance of the adjusted points should be smaller than that of the original points. So that, Eq.
If the B-spline curve fails to meet the precision requirement after several iterations, the correction process which includes knot insertion and foot point parameter (FPP) upgrading should be involved [23]. To accelerate the iteration, new knots and corresponding control points are inserted into intervals where the fitting tolerance is exceeded. FPP upgrading means to replace the tool tip point parameters with the FPPs of tool tip points so that the distance vector k can be estimated more precisely. In this paper, knot insertion and FPP upgrading are conducted every 10 iterations unless the terminating condition has been satisfied.
The knot insertion is mainly composed of two steps: evaluating the degrees of freedom (DoFs) of knot intervals and searching the insert interval. The DoF of each control point in the kth iteration can be evaluated by t j � ‖ exceeds the fitting tolerance, the corresponding knot interval and its DoF will be found and evaluated. The DoF of the knot interval u a , u a+1 in the kth iteration can be evaluated by The interval with the maximum k a will be chosen as the insert interval, e.g., u ins a , u ins a+1 , and the insert knot is u ins = u ins a +u ins a+1 2 . The FPP of each tool tip point j can be calculated by Eq. (15). However, there may exist several foot points on the curve near one tool tip point. In the case shown in Fig. 7, the result of Newton's iteration may converge to an inappropriate value, which can result in the deterioration of fitting precision. Thus, the identification and correction of inappropriate FPPs should be conducted to improve the stability of the algorithm.
Since the initial FPPs are the tool tip point parameters t j n j=0 , the FPP sequence t foot j n j=0 is monotonically increasing. Therefore, an FPP t foot j is considered as inappropriate when the following condition is satisfied.
Note that the initial value can greatly influence the result of Newton's iteration. Then, for an inappropriate result, the FPP is calculated repeatedly with different initial values from the sequence t j−w , ⋯ , t j−1 , t j+1 , t j+w , which usually t j−1 , t j+1 is enough. After that, the parameter of the nearest foot point to   In addition, if a closed-fitted B-spline curve is desired, endpoint constraints should be considered. To ensure that the fitted B-spline curve interpolates the endpoints, control points 0 and m will not be moved in each iterative step, i.e., 0 = 0 and m = n . To guarantee the G1 continuity at the endpoints, control points 1 and m−1 will be predetermined as follows and unchanged during the iteration.
where V is the unit tangent vector at the endpoints.
During the chord error refinement of ELSPIA, the control point positions of the fitted B-spline curve are optimized to minimize the objective function in Eq. (22). The iterative steps of ELSPIA during chord error refinement can be described as where can be estimated through = 2 min + max with the largest and the smallest eigenvalues of T + , and is a diagonal matrix defined as In each iterative step, a knot interval u a , u a+1 where the chord error exceeds the tolerance is first to be determined. Secondly, a new knot will be inserted in the interval, and a new control point will be added meanwhile. Then three control points a−2 , a−1 , and a corresponding to u a , u a+1 will be adjusted, and other control points remain unchanged. Finally, the iteration stops when the chord errors fulfill the tolerance requirement. The method of chord error estimation will be discussed in the next section.

Approximation error estimation
The approximation error estimation is composed of fitting error estimation and chord error estimation. The fitting errors denote the deviations between the original tool tip points and the fitted B-spline curve, which can be estimated by calculating the distance from each tool tip point to its foot point through Eqs. (15) and (16). The   chord errors indicate the maximum deviations between linear segments and the fitted B-spline curve. A Newton iteration method based on geometric relationships to compute the chord errors is proposed in Ref. [23].
To shorten the computing time, the chord errors are estimated as follows in this paper. Firstly, the fitted B-spline curve is divided into curve segments by the foot points of the original tool tip points. Secondly, for each curve segment, a group of points is sampled with uniformly distributed parameters. Finally, calculate the distances between the sample points and the corresponding linear segments. For each group of sample points, the maximum distance will be regarded as the chord error between the corresponding linear segment and the B-spline curve. The schematic diagram of chord error estimation is shown in Fig. 8.

Experimental study
In this section, curve fitting experiments are conducted to prove the effectiveness of the proposed toolpath smoothing method on maximum curvature reduction. The proposed method generates smooth B-spline toolpaths based on the ELSPIA method with tool tip position adjustment. The contrast method is the existing ELSPIA method without tool tip position adjustment. Furthermore, interpolation experiments are conducted on a commercial CNC system to validate the effectiveness of the proposed method in improving machining efficiency.

Curve fitting experiment results
Curve fitting experiments, performed on the butterfly and rabbit toolpaths, are implemented off-line by MATLAB 2018a on a personal computer with a 64-bit operating system, a RAM of 16 GB, and a CPU of AMD Ryzen 7 4800H 2.9 GHz. The butterfly toolpath and the rabbit toolpath consist of 702 and 836 tool tip points, respectively. The parameters of tool tip position adjustment are listed in Table 1. Consequently, 112 out of 702 tool tip points of the butterfly toolpath are adjusted with a computation time of 1.97 s. For the rabbit toolpath, 161 out of 836 tool tip points need adjustment, and the computation time is 5.26 s.
The fitting tolerance tol is set to be 0.03 mm, and the neighborhood parameter s of dominant point selection is 3 for the butterfly toolpath and 4 for the rabbit toolpath. As a result, 215 and 206 dominant points are selected from the original butterfly toolpath and rabbit toolpath. After tool tip position adjustment, 215 and 208 dominant points are chosen from the adjusted toolpaths, which demonstrates that the proposed method preserves the overall shape of the toolpaths well.
The fitted B-spline curves of the butterfly and rabbit toolpaths are shown in Figs. 9 and 10. A qualitative illustration of the fitted butterfly curve geometry is provided. In flat areas of the butterfly toolpath and the rabbit toolpath, detailedly displayed in Figs. 9b and 11b, fitted curves obtained by different methods are close to each other with little deviation. However, in a region near a sharp corner of the butterfly toolpath, the resultant curve of the proposed method tends to be less curved than that of the contrast method, which is demonstrated in Fig. 9c to f. For the rabbit toolpath, the B-spline curve generated by the proposed method also turns out to be more flat near sharp corners, which is demonstrated in Fig. 11c to f. The approximation errors of the butterfly and rabbit toolpaths are shown in Figs. 10 and 12. It can be observed that both fitting errors and chord errors of the B-spline curves generated by the proposed and contrast method are within the fitting tolerance.   Figs. 13b and 14b that the curvature peak values of the curves generated by the proposed method are much less than those of the contrast method. Furthermore, the local curvature maximums at eight different positions are selected to be compared for both cases, which are listed in Table 2. In the case of the butterfly curve, seven out of the eight local curvature maximums of the fitted curve produced by the proposed method are less than those of the contrast method. The maximum reduction of the local curvature maximums is 59.13%, from 32.71 to 13.37. In the case of the rabbit curve, all of the eight local curvature maximums obtained by the proposed method are less than those of the contrast method, with a maximum reduction of 42.86%, from 1.19 to 0.68.
To sum up, it has been validated qualitatively and quantitatively that fitting toolpaths with the proposed tool tip

Interpolation experiment results
Feedrate scheduling and interpolation experiments are performed on the 5-axis platform shown in Fig. 15. The platform is composed of a SINUMERIK 840D sl CNC system, SIMOTIN D-425 motion control system, SINAM-ICS S120 drive system, and servo motors. The parameters of feedrate scheduling and interpolation are listed in Table 3.
The feedrate profiles after feedrate scheduling and interpolation are shown in Fig. 16. Compared with the contrast method, the machining time of the butterfly curve reduces from 6.830 to 6.652 s, and the machining efficiency is increased by 2.61% with the proposed method. In addition, the machining efficiency of the rabbit curve is increased by 2.37%, with a reduction of the machining time from 7.944 to 7.756 s. As a result, the effectiveness of the proposed method in improving machining efficiency has been verified.
The deviations of the interpolated points, output by the CNC system, from the original toolpaths are calculated and demonstrated in Fig. 17. As a result, the deviation of the proposed method is at the same level as that of the contrast method for both the butterfly and rabbit toolpaths. In addition, the deviations do not exceed the fitting tolerance of 0.03 mm. It proves that the proposed method maintains machining accuracy while improving efficiency.

Conclusion
B-spline curves are widely used for linear toolpath smoothing to obtain smooth feedrate profiles and high machining efficiency. Most research on B-spline curve approximation lay much emphasis on constraining the approximation errors. Different from most existing methods, a global toolpath smoothing method based on B-spline approximation with tool tip position adjustment which aims at maximum curvature reduction is proposed in this paper. To decrease the maximum curvatures of fitted curves, toolpaths are first split into several subdivisions, and the subdivisions that need to be adjusted are identified according to the discrete curvatures of tool tip points. For each subdivision that needs adjustment, tool tip positions are optimized to minimize the maximum discrete curvature. The new toolpaths after adjustment are fitted by an improved ELSPIA method. The proposed approach has the following advantages: (a) the overall shape of the toolpath is well-preserved after tool tip position adjustment.
(b) Global G2 continuity of the fitted B-spline curve is realized. (c) The maximum curvature of the fitted curve is significantly decreased, and the machining efficiency is improved. (d) Approximation errors are strictly confined to tolerance. (e) The tool tip position adjustment can be flexibly combined with not only the ELSPIA, but also other different approximation methods. Although it is proposed for 3-axis machine tools, the method can be extended to apply to 5-axis machine tools. The experimental results demonstrate that the proposed method can decrease the maximum curvature by about 50% and, as a result, improve the machining efficiency by about 3%.
Author contribution All authors contributed to the study's conception and design. Material preparation, data collection, and analysis were performed by Li Hua, Nuodi Huang, and Bowen Yi. The first draft of the manuscript was written by Li Hua, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.