A Trajectory Planning Method for Polishing Optical Elements Based on a Non-Uniform Rational B-Spline Curve

Abstract

Computer-controlled optical surfacing (CCOS) has been successfully applied in industrial production, as it can precisely correct the surface error by converting the dwell time of the polishing tool into the feed-rate along the polishing trajectory. [...]a trajectory planning method is the key factor affecting high-quality, high-efficiency polishing. [...]the SOR iterative algorithm [28] is used to find stable numerical solutions of Equation (5) as described below. According to affine invariant principle of NURBS [31], the four-dimensional (4D) space constructed by the control points and weight factors can be expressed as: Cω(u)=∑i=1n+1Ni,3(u)[ωi diωi]=∑i=1n+1Ni,3(u)diω Then, the NURBS defined by Equation (1) can be regarded as the projection of the curve Cω(u) in 4D space on the centre of the hyperplane ω=1 . Step 4: the parameter interval (0,1) was sectioned by the equivalent distance up−low to generate the knot vector (u0,u1,…,un) .

1. Introduction

With the rapid development of astronomy, space exploration, and advanced optical instruments, optical elements are being increasingly widely used. The application demands for high-quality and high-efficiency elements present distinct higher requirements for the process technology used in such elements [1]. Computer-controlled optical surfacing (CCOS) has been successfully applied in industrial production, as it can precisely correct the surface error by converting the dwell time of the polishing tool into the feed-rate along the polishing trajectory. Therefore, a trajectory planning method is the key factor affecting high-quality, high-efficiency polishing.

Despite the extent of research on path planning in some fields [2,3,4], investigation of optical polishing has been rather limited. Although the problem of discontinuous surfaces between the adjacent mm-sized short line segments has attracted wide concern when using parametric curve theory [5,6,7], frequent acceleration and deceleration result in poor realisation precision of dwell time, i.e., the runtime of trajectory (hereafter referred to as runtime), which further influences the polishing quality of elements and the convergence rate of surfaces [8]. The main focus of this paper is to investigate an interpolation scheme to overcome the above problem, which has two main progressive aspects: fitting and interpolation of parametric curves.

The fitting of parametric curves is a process of converting discrete short line segments into parametric curves. At present, many scholars have carried out research based on the dominant points using a non-uniform rational B-spline (NURBS). Park [9,10] proposed a method for determining the dominant points according to the discrete curvature. Zhou [11] and Xu [12] improved this method by taking concave-convex turning points and extreme points on the curvature curve as dominant points. To improve the fitting precision, Zhao [13] proposed curve fitting taking squared distance minimisation (SDM) as the evaluation index. Although the dominant points based method is easy with regard to calculation and interpolation, it may result in the loss of runtime at non-dominant points. Thus, only the global fitting method is suitable to the optical polishing. To do so, Yang [14] proposed an optimisation algorithm by establishing the evaluation function for deviation of the fitted distance. Li [15] and Lin [16] classified the trajectory into the different forms of NURBS and then employed a piecewise fitting method for real-time implementation. Based on Gaussian elimination and the continuous short block (CSB) look-ahead algorithm, Tsai [17] and Wang [18] realised the on-line transformation from short line segments to NURBS: however, for the global fitting method, the simplified strategy that setting all the weight factors as 1 eliminates the regulating effect on the fairness of curves and easily causes curvature saltation on the trajectory. Therefore, generating a fair trajectory based on NURBS is the first problem facing the polishing trajectory planning technique.

The interpolation of parametric curves is a process that discretises the NURBS into numerical control (NC) commands. Speed planning, as the critical step in the interpolation process, has been an area of research for numerous scholars: this can be classified into the time-optimal approach and the non-time-optimal approach. The time-optimal approach deals with the planning problem by taking the minimisation of the motion time as the objective to promote manufacturing efficiency [19,20]. For example, Timar [21] proposed a speed planning scheme for NC interpolation by the use of the optimal control theory. On this basis, Sencer [22] and Lu [23] considered the driving capacity constraint and trajectory precision constraint in the interpolation, respectively. The non-time-optimal approach usually deals with the planning problem by taking the minimisation of speed fluctuations as the objective [24]. Various methods, such as the feed-rate evolutionary algorithm [25], the equidistance quaternion method [26], and the improved Adams-Malton algorithm [27] are employed to decrease speed fluctuations and guarantee steady, continuous-trajectory operation. Although the effectiveness of the method has been validated experimentally, is cannot be applied directly to optical polishing because the effect of acceleration and deceleration on the realisation precision of runtime is not considered.

Driven by the practical needs to improve the quality of optical polishing, this paper presents a systematic trajectory planning method that particularly enhances the realisation precision of runtime. Following this introduction, Section 2 calculates all the control points of NURBS through numerical solution to overcome the problem of calculation efficiency. In Section 3, an optimisation algorithm of fairing NURBS is established by taking the shear jerk of trajectory as an evaluation index. Section 4 then proposes an interpolation scheme with which to minimise the realisation error of runtime by planning the feed-rate of trajectory according to the given runtime between adjacent NC codes. Section 5 reports experiments on a prototype machine which shows that the proposed trajectory planning method is more accurate than the linear interpolation method. Conclusions are drawn in Section 6.

2. Trajectory Fitting Based on the NURBS Curve

In this section, according to the given NC codes, all the control points of NURBS are solved, which provides the necessary mathematical model for the interpolation of parametric curves. The basic settings are as follows:

(1) Considering the stability, ease of use and calculation efficiency, cubic NURBS is employed as the fitting tool. (2) The uniform parametric method is used for trajectory fitting, because the polishing elements have a large radius of curvature and the chord lengths between NC codes are distributed uniformly. (3) All weight factors are set to 1. In this case, the NURBS can be treated as a cubic B-spline to simplify the calculation process.

According to the basic theory of the NURBS, a segment of NURBS C(u)(0≤u≤1) can be determined based on the four adjacent control points _dk(k=0,1,2,3) . As shown in Figure 1, C(0) and C(1) separately refer to the start point and the end point of the NURBS segment. Based on the aforementioned setting, the NURBS segment can be directly written as:

C(u)=16[^u3u2u1][−13−313−630−30301410][_d0d1d2d3]

Thus, point _Ci can be expressed as:

_Ci=16(_di−1+4_di+_di+1)

where _Ci denotes the ith NC code in the n+1 lines of NC codes and also is taken as the start point of the ith NURBS segment. _di−1,_di,_di+1 respectively denote the three control points corresponding to _Ci . To calculate _C1 and _Cn+1 , it is defined that on the condition that i<1 , then _di=_d1 and on the condition that i>n+1 , then _di=_dn+1 . Equations (3) and (4) can thus be obtained:

_C1=16(_d1+4_d1+_d2)=56_d1+16_d2

_Cn+1=16(_dn+4_dn+1+_dn+1)=16_dn+56_dn+1

Rewriting Equations (2)–(4) in matrix notation yields:

Ad=CA=[51000…014100…001410…0⋮⋮⋮⋮⋮…⋮00001410000015], d=[_d1d2d3⋮_dndn+1],C=[6_C16_C26_C3⋮6_Cn6_Cn+1]

If the traditional solution methods such as Gauss elimination method and LU decomposition are directly applied to the Equation (5), it may lead to some undesirable phenomena (such as excessive calculation time) due to the large quantity of NC codes for polishing. Hence, the SOR iterative algorithm [28] is used to find stable numerical solutions of Equation (5) as described below.

The coefficient matrix A can be divided into:

A=D−L−U

where A=diag(_a11,_a22,…,_an+1n+1) , −L denotes a strictly lower triangular matrix, whose elements below the principal diagonal are corresponding elements of A . −U denotes a strictly upper triangular matrix, whose elements above the principal diagonal are corresponding elements of A . Owing to the matrix D being invertible, Equation (5) can be modified to:

d=^D−1(L+U)d+^D−1C

In this case, corner marks (k) and (k+1) are added in Equation (7) to identify the number of iterations. Then, the elementary iterative scheme of d can be written as:

^d(k+1)=_M1 ^d(k)+^D−1C

where _M1=^D−1(L+U) . The component form of Equation (8) can be expressed as:

_di(k+1)=1_aii(−∑j=1j≠in+1_{aij dj(k)}+_Ci), i=1,2,…,n+1

which is also known as the Jacobi iteration scheme, noting that, before calculating _di(k+1) , the iterative values of the first i−1 components in ^d(k+1) have been generated, which are more approximate to the true value than the results obtained by the previous iteration. Therefore, _d1(k),_d2(k),…,_di−1(k) can be replaced by _d1(k+1),_d2(k+1),…,_di−1(k+1) to make _di(k+1) closer to the true value. To improve the convergence rate of the iteration further, a proper parameter μ is selected to conduct the weighted averaging on the aforementioned iterative scheme:

_di(k+1)=μ_d˜i(k+1)+(1−μ)_di(k),i=1,2,…n+1

with:

_d˜i(k+1)=1_aii(−∑j=1i−1_{aij dj(k+1)}−∑j=i+1n+1_{aij dj(k)}+_Ci)

Substituting Equation (11) into Equation (10) leads to the SOR iteration scheme:

_{aii di(k+1)}+μ∑j=1i−1_{aij dj(k+1)}=_{aii di(k)}−μ∑j=in+1_{aij dj(k)}+μ_Ci,i=1,2,…,n+1

Note that A is a tridiagonal positive definite matrix, so the optimal value _μopt of the relaxation factor can be expressed as [29]:

_μopt=21+1−^[ρ(_M1)]2

where ρ(_M1) denotes the spectral radius of _M1 .

It is worth noting that, for a flat element, the aforementioned method can be directly applied to calculate the NURBS trajectory, but for a curved element, the polishing shaft should always lie along the normal direction of the element. In this case, it is necessary to solve two trajectories of end-point and reference-point on the polishing shaft to determine the polishing attitude (for more details, please see [21]).

3. Fairing of the NURBS

Although the constructed cubic NURBS by the above method satisfies the G2 continuity characteristics, i.e., the second-order derivative functions of the trajectory is continuous, the motion stability of the trajectory is still influenced by curvature saltation. In this section, a fairing optimisation method is proposed by adjusting the weight factors of NURBS.

The fairing optimisation can be classified into two categories, i.e., global and local fairing according to the number of the adjusted point on NURBS. Considering the significant computational burden of global fairing caused by the large quantity of NC codes for polishing, it is more reasonable to modify the outlier points with curvature saltation by the use of local fairing, which are selected from all NC codes [30]. A filtering process is needed to eliminate the influence of curvature fluctuations caused by discrete calculation.

In the fairing optimisation, the shear jerk of the outlier point is taken as the evaluation index:

_Dj=_κnext1−_κj‖_Cnext1−_Cj‖−_κj−_κlast1‖_Cj−_Clast1‖

where _κj denotes the curvature at the jth outlier point _Cj , _κnext1 and _κlast1 denote the curvature of the _Cnext1 and _Clast1 , which are on two adjacent sides of the outlier point. The index indicates the curvature changes of the adjacent outlier points.

To generate the weight factors of the various outlier points, the objective function is defined as:

L=∑j=1k_Dj2+^{(_Cj−_Cj,0)2}

where _Cj,0 denotes the jth outlier point before optimisation. It can be seen from Equation (15) that the objective function is composed of two parts: part one is the curvature changes of the outlier point after optimisation, and part two is the adjustment amplitude of outlier points before, and after, optimisation. Thus, the objective function means that fairing optimisation is performed on the premise of modifying the NURBS as little as possible.

According to affine invariant principle of NURBS [31], the four-dimensional (4D) space constructed by the control points and weight factors can be expressed as:

^Cω(u)=∑i=1n+1_Ni,3(u)[_{ωi diωi}]=∑i=1n+1_Ni,3(u)_diω

Then, the NURBS defined by Equation (1) can be regarded as the projection of the curve ^Cω(u) in 4D space on the centre of the hyperplane ω=1 . Based on Equation (2), it can be seen that:

_Cjω=16_dj,last1ω+23_djω+16_dj,next1ω

_Cj,last1ω=16_dj,last2ω+23_dj,last1ω+16_djω

_Cj,next1ω=16_djω+23_dnext1ω+16_dnext2ω

where _djω denotes the control point corresponding to the jth outlier point and _dnext1ω,_dnext2ω and _dlast1ω,_dlast2ω denote the control points on the two adjacent sides of the control point _djω , respectively.

Furthermore, Equation (15) can be written in 4D space as:

^Lω=∑j=1k^(_Djω)2+^{(_Cjω−_Cj,0ω)2}

where _Cjω−_Cj,0ω can be equivalently simplified as _djω−_dj,0ω and _Djω can be approximately expressed as [32]:

_Djω≈^{(_Cj,next1ω)′​′}−^{(_Cjω)′​′}_lj,next1ω−^{(_Cjω)′​′}−^{(_Cj,last1ω)′​′}_lj,last1ω

where ^{(_Cjω)′​′} denotes the second derivative of the NURBS at _Cjω , _lj,next1ω=‖_Cj,next1,0ω−_Cj,0ω‖ and _lj,last1ω=‖_Cj,0ω−_Cj,last1,0ω‖ . Substituting Equation (21) into Equation (20) yields:

^Lω=∑j=1k(^{(_Cj,next1ω)′​′}−^{(_Cjω)′​′}_lj,next1ω−^{(_Cjω)′​′}−^{(_Cj,last1ω)′​′}_lj,last1ω)+^{(_djω−_dj,0ω)2}

To calculate the minimum value of ^Lω , the partial derivative of Equation (22) about _djω is set to 0:

∂^Lω∂_djω=∑j=1k2(_dj,next2ω−3_dj,next1ω+3_djω−_{dj,last1ωlj,next1ω}−_dj,next1ω−3_djω+3_dj,last1ω−_{dj,last2ωlj,last1ω})·(3_lj,next1ω+3_lj,last1ω)+2(_djω−_dj,0ω)=0

Then, the equations for _d1ω,_d2ω,…,_dkω can be expressed as:

^{Aω dω}=^CωAω=diag(^{(3_lj,next1ω+3_lj,last1ω)2}+1)^dω=^{(_d1ω,_d2ω,…,_dkω)TCω}=^{(_C1ω,_C2ω,…,_Ckω)T}_Cjω=_dj,0ω−(_dj,next2ω−3_dj,next1ω−_{dj,last1ωlj,next1ω}−_dj,next1ω+3_dj,last1ω−_{dj,last2ωlj,last1ω})(3_lj,next1ω+3_lj,last1ω)

According to Equation (24), the optimised control points and weight factors corresponding to the outlier points can then be generated.

4. NURBS Interpolation

The NURBS interpolation is used to discretise the parametric curve to the NC commands based on the planned feed-rate. In this section, an interpolation method for optical polishing is proposed aiming to minimise the realisation error of the trajectory runtime.

4.1. Feed-Rate Planning

Feed-rate planning is the main influencing factor in interpolating the NC commands along the trajectory. For the optical polishing, to guarantee the desired runtime of trajectory, the specific method is displayed as follows:

Step 1: considering that now there is no analytical solution to calculate the length of NURBS, the Simpson formula is used to obtain the estimation of the length through numerical iteration:

p=up−low6(f(low)+4f(mid)+f(up))

with:

f(u)=dsdu=^x′(u)+^y′(u)+^z′(u)

where ^x′(u),^y′(u),^z′(u) respectively denote the first-order derivatives of x(u),y(u),z(u) , which are the one-dimensional curves of C(u) along x,y,z axes. up and low denote the upper and lower boundaries of u , mid=(up+low)/2 .

Step 2: the length between the two points corresponding to the parameters up and low is calculated:

l=‖C(up)−C(low)‖

Step 3: the error between the aforementioned two lengths is calculated:

e=|p−l|

The convergence threshold [e] is given. If e>[e] , let up=mid and repeat Steps 1 and 2. If e<[e] , turn to Step 4.

Step 4: the parameter interval (0,1) was sectioned by the equivalent distance up−low to generate the knot vector (_u0,_u1,…,_un) . Thus, the length of the NURBS is:

s≈∑i=1n‖C(_ui)−C(_ui−1)‖

Equipped with the length at hand, the S-curve motion law is invoked to guarantee that the feed-rates of the adjacent NURBS segments are changed smoothly. Moreover, the initial sections of the trajectory segments are defined as the feed-rate transition zone. Then, the feed-rates remain constant until the end of the trajectory segments. In this case [33]:

s=2_{vlow ta}+J_ta3+_{vup tv}

_vup=_vlow+J_ta2

where _vlow and _vup denote the initial and final feed-rates of the trajectory segment. J , 2_ta and _tv denote the jerk, acceleration (deceleration) time and the uniform motion time, respectively. Thus, the runtime of the trajectory segment is _td=2_ta+_tv .

It is worth noting that previous studies have shown that the feed-rates, when limited by the runtime of the trajectory segments, are much lower than the maximum value which is constrained by the chord error and the driving capacity. Therefore, there is no need to check the feed-rate again.

4.2. NURBS Interpolation

The essence of interpolation is to generate the NC command along the trajectory according to the period _ts . As each interpolation point of the NURBS corresponds to one curve parameter, only the curve parameters need to be solved.

Taking u as the function of t , the second-order Taylor expansion can be expressed as:

_ui+1=_ui+_{dudt|t=ti ts}+12_{^d2ud^t2|t=ti ts2}

The feed-rate v(u) can be written as:

v(u)=‖dC(u)du‖=‖dC(u)dududt‖

Furthermore:

dudt=v(u)‖^C′(u)‖

Calculating the derivative of Equation (34):

^d2ud^t2=dv(u)dududt‖^C′(u)‖−v(u)d(‖^C′(u)‖)dududt^{‖C′(u)‖2}

where:

d(‖^C′(u)‖)du=^C″(u)⋅^C′(u)‖^C′(u)‖

Substituting Equations (33)–(35) into Equation (32) gives:

_ui+1=_ui+_vC(_ui)‖^C′(_ui)‖_ts+12[ _^v′C(_ui)^{‖C′(_ui)‖2}_vC(_ui)−^C″(u)⋅^C′(u)^{‖C′(u)‖4}_vC2(_ui)]_ts2

Substituting _ui+1 into the curve equation obtained through the fairing optimisation described in Section 3, the next interpolation point can be acquired. Repeating this process until:

_ti=ceil(_tdts)*_ts

where ceil(·) denotes an integer that is rounded up. In this way, interpolation of all trajectory segment can be completed.

5. Experiments

Both simulation and experiments were carried out to validate the effectiveness of the presented method on the prototype of the hybrid polishing robot. As shown in Figure 2, it is mainly composed of a 6-DOF (degrees of freedom) hybrid robot, a polishing effector, a magnetic worktable, a column, and a CNC system. The hybrid robot is composed of a 3-DOF (3UPS and UP) parallel mechanism and a 3-DOF wrist. The UP limb and the wrist form a UPS or UPRRR limb. Here, R, U, S, and P represent, respectively, revolute, universal, spherical, and prismatic joints, and the underlined P, S, and R denote the actuated prismatic, spherical, and revolute joints, respectively. The CNC system is built upon an IPC+PMAC open architecture, consisting of a host control computer responsible for reconstruction of the parametric curves, trajectory interpolation, and NC command generation, and a PMAC motion controller for servo-control of the actuated joints.

Without loss of generality, a segment of NC codes was taken from the polishing trajectory to validate the effectiveness of the proposed interpolation method. According to the SOR iterative scheme mentioned in Section 2, Figure 3 shows the two fitted trajectories. Then, the optimisation algorithm described in Section 3 is used to smooth the curvature of the trajectory. The change threshold of curvature for judging outlier points is given as 0.01. As shown in Table 1, the fairness of the two trajectories is both significantly improved. The maximum curvatures of the two trajectories are reduced by as much as 79.7% and 63.3% and the maximum absolute values of shear jerks are decreased by 91.2% and 90.2%, correspondingly. Considering the optimisation and interpolation methods for the two trajectories are same, experimental results of the end-point trajectory are just shown in the following discussion for the sake of simplicity.

Based on the optimisation trajectory shown in Figure 4, the discrete NC command sequences were generated and sent to the PMAC motion controller, which were mapped into the servo-command of actuated joints through an inverse kinematic model. The PMAC motion controller synchronously gathered the positions and velocities fed back from the servo-motors. The interpolation and sampling periods are 10 and 20 ms, respectively. To validate the effectiveness of the method, a comparison experiment was carried out utilizing the linear interpolation method. Figure 5 shows the experimental result, which is computed by the feedback positions of all the actuated joints. It can be seen from the partial enlarged view that the NC commands generated by the proposed method are closer to the NC codes than those found using linear interpolation, which means that the removal position of polishing process can be reached more accurately. Figure 6 and Figure 7 show the runtime between the adjacent NC nodes. Compared with the desired value, the proposed method is able to realise the runtime more precisely, which indicates that the removal quantity during the polishing process can be more precisely controlled, correspondingly. To evaluate the effect of the interpolation method, the indices are defined as the interpolation error ( _ei ) and the runtime error of the trajectory ( _et ):

_ei=_{‖Ci0−Ci‖2},i=1,2,…,n

_eti=|_ti0−_ti|,i=1,2,…,n

where _Ci0 and _Ci denote the ith desired NC code and actual NC code after interpolation, _ti0 and _ti denote the desired runtime and the actual runtime between the i − 1th and ith NC codes. It can be seen from Table 2 that, compared to the linear interpolation method, the interpolation error and the runtime error generated through the proposed method are an order of magnitude smaller. These data further imply that the proposed interpolation method has more precision than the linear interpolation method in the polishing process.

To validate the modification effect of the interpolation method on the surface error of optical elements, polishing experiments were carried out on fused silica elements by, respectively, using the linear interpolation method and the proposed interpolation method, as shown in Figure 8. Using the Nanovea contour graph, the polished zone (70 mm × 70 mm) of the elements was detected. The experimental results obtained using the two interpolation methods are displayed in Figure 9 and Figure 10, respectively. It can be seen, from the figures, that, after conducting linear interpolation-based polishing, the surface error (PV) decreases from 11.894λ(λ=633nm) to 7.422λ . In contrast, using the proposed interpolation method, the surface error (PV) is reduced from 11.282λ to 6.267λ . The convergence rate of the surface error is increased from 37.59% to 44.44%, which further verifies the effectiveness of the proposed method.

6. Conclusions

A new trajectory planning method for optical polishing is proposed in this paper. It is developed to generate NC commands that can decrease the interpolation error and the runtime error. First, to obtain the NURBS trajectory without loss of the information about the runtime, a global fitting method is presented based on SOR iteration theory to deal with the problem of the computational burden arising from the large quantity of NC codes for polishing. Then, taking the shear jerk of the NURBS as the evaluation index, a fairing optimisation method is carried out to smooth the curvature saltation of the trajectory. Finally, the feed-rate planning method and the path interpolation scheme are proposed to reduce the realisation error of the trajectory runtime.

Simulation results verify that the fairing optimisation proposed in this research can modify the curvature saltation at the expense of trajectory accuracy compared to the given NC codes; however, the curvature saltation only occurs at the turning point of the trajectory and the trajectory error magnitude generated by the optimisation algorithm is consistent with the linear interpolation. Therefore, the trajectory accuracy is not treated as the index with which to evaluate the smoothing effect in Section 5. The effect of the optimisation algorithm on the optical polishing will be investigated in future work.

It can be seen from Figure 5, Figure 6 and Figure 7, and Table 2, that the runtime error of the proposed interpolation method arises because the trajectory runtime cannot be divided exactly by the interpolation period and is rounded up to an integer. In contrast, the runtime error of the linear interpolation method arises as a result of the acceleration and deceleration that occurs frequently between adjacent short line segments. Similar to the runtime error, the interpolation error of the proposed interpolation method is a result of the numerical calculation error of the trajectory length and that of the linear interpolation method is due to the transition for the discontinuity of the adjacent short line segments. Then the convergence rate of the surface error is increased with the help of the improvement in the interpolation, and runtime, errors, validating the effectiveness of the proposed interpolation process. For the aforementioned error, their individual effects on the optical polishing need to be further indicated in future work.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

Trajectory	Before		After
	Maximum Curvature	Maximum Shear Jerk	Maximum Curvature	Maximum Shear Jerk
End-point	1.249	0.285	0.253	0.025
Reference-point	1.380	0.325	0.506	0.032

Method	Interpolation Error (mm)		Runtime Error (s)
	Maximum	Mean	Maximum	Mean
Proposed	0.005	0.002	0.010	0.005
Linear	0.091	0.076	0.130	0.076

Author Contributions

D.Z. wrote the manuscript and performed the experiments, H.G. analysed the data and revised the manuscript. All authors discussed the results and commented on the manuscript.

Funding

This research was partially funded by the National Science and Technology Major Project of China (grant no. 2017ZX04022001-206), the National Natural Science Foundation of China (grant no. 51420105007), and EU H2020-MSCA-RISE 2016 (grant no. 734272).

Conflicts of Interest

The authors declare no conflict of interest.

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Rupp, W. Conventional optical polishing techniques. Int. J. Opt. 1971, 18, 1-16.

2. Kaltsoukalas, K.; Makris, S.; Chryssolouris, G. On generating the motion of industrial robot manipulators. Rob. Comput. Integr. Manuf. 2015, 32, 65-71.

3. Makris, S.; Tsarouchi, P.; Matthaiakis, A. Dual arm robot in cooperation with humans for flexible assembly. CIRP Ann. 2017, 66, 13-16.

4. Lei, W.; Wang, S. Robust real-time NURBS path interpolators. Int. J. Mach. Tools Manuf. 2009, 49, 625-633.

5. Jahanpour, J.; Alizadeh, M. A novel acc-jerk-limited NURBS interpolation enhanced with an optimized S-shaped quintic feedrate scheduling scheme. Int. J. Adv. Manuf. Technol. 2015, 77, 1889-1905.

6. Liu, X.; Peng, J.; Si, L. A novel approach for NURBS interpolation through the integration of acc-jerk-continuous-based control method and look-ahead algorithm. Int. J. Adv. Manuf. Technol. 2017, 88, 961-969.

7. Zhao, J.; Li, L.; Wang, G. Research of NURBS Curve Real-time Interpolation Feed Speed Planning. Tool Eng. 2015, 49, 7-11.

8. Lin, M.; Tsai, M.; Yau, H. Development of a dynamics-based NURBS interpolator with real-time look-ahead algorithm. Int. J. Mach. Tools Manuf. 2007, 47, 2246-2262.

9. Park, H.; Lee, J. B-spline curve fitting based on adaptive curve refinement using dominant points. Comput.-Aided Des. 2007, 39, 439-451.

10. Park, H. B-spline surface fitting based on adaptive knot placement using dominant columns. Comput.-Aided Des. 2011, 43, 258-264.

11. Zhou, H.; Wang, Y.; Liu, Z. On non-uniform rational B-splines curve fitting based on the least control points. J. Xian Jiaotong Univ. 2008, 42, 73-77.

12. Xu, J. B-spline Curve Approximation Based on Feature Points Automatic Recognition. J. Mech. Eng. 2009, 45, 212-217.

13. Zhao, H.; Lu, Y.; Zhu, L. Look-ahead interpolation of short line segments using B-spline curve fitting of dominant points. P. I. Mech. Eng. B-J. Eng. 2014, 229, 1-13.

14. Yang, X.; Hu, Z.; Zhong, Z. Research on the NURBS Curve Fitting for Tool Path Generation. China Mech. Eng. 2009, 20, 983-984.

15. Li, W.; Liu, Y.; Yamazaki, K. The design of a NURBS pre-interpolator for five-axis machining. Int. J. Adv. Manuf. Technol. 2008, 36, 927-935.

16. Lin, K.; Ueng, W.; Lai, J. CNC codes conversion from linear and circular paths to NURBS curves. Int. J. Adv. Manuf. Technol. 2008, 39, 760-773.

17. Tsai, M.; Nien, H. Development of a real-time look-ahead interpolation methodology with spline-fitting technique for high-speed machining. Int. J. Adv. Manuf. Technol. 2010, 47, 621-638.

18. Wang, J.; Yau, H. Real-time NURBS interpolator: application to short linear segments. Int. J. Adv. Manuf. Technol. 2009, 41, 1169-1185.

19. Timar, S.; Farouki, R. Time-optimal traversal of curved paths by Cartesian CNC machines under both constant and speed-dependent axis acceleration bounds. Rob. Comput. Integr. Manuf. 2007, 23, 563-579.

20. Sencer, B. Smooth Trajectory Generation and Precision Control of 5-Axis CNC Machine Tools. Ph.D. Thesis, The University of British Columbia, Vancouver, BC, Canada, October 2009.

21. Boyadjieff, C.; Farouki, R.; Timar, S. Smoothing of Time-Optimal Feed rates for Cartesian CNC Machines; Springer: Berlin, Germany, 2005; pp. 84-101.

22. Sencer, B.; Altintas, Y.; Croft, E. Feed optimization for five-axis CNC machine tools with drive constraints. Int. J. Mach. Tools Manuf. 2008, 48, 733-745.

23. Lu, L.; Zhang, L.; Ji, S. An offline predictive feedrate scheduling method for parametric interpolation considering the constraints in trajectory and drive systems. Int. J. Adv. Manuf. Technol. 2016, 83, 2143-2157.

24. Zhao, H.; Zhu, L.; Ding, H. A parametric interpolator with minimal feed fluctuation for CNC machine tools using arc-length compensation and feedback correction. Int. J. Mach. Tools Manuf. 2013, 75, 1-8.

25. Sun, Y.; Zhao, Y.; Bao, Y. A novel adaptive-feedrate interpolation method for NURBS tool path with drive constraints. Int. J. Mach. Tools Manuf. 2013, 77, 74-81.

26. Zhang, J.; Zhang, L.; Zhang, K. Double NURBS trajectory generation and synchronous interpolation for five-axis machining based on dual quaternion algorithm. Int. J. Adv. Manuf. Technol. 2016, 83, 2015-2025.

27. Peng, J.; Liu, X.; Si, L. A Novel Approach for NURBS Interpolation with Minimal Feed Rate Fluctuation Based on Improved Adams-Moulton Method. Math. Probl. Eng. 2017, 2017, 1-10.

28. Bai, Z.; Parlett, B.; Wang, Z. On generalized successive over relaxation methods for augmented linear systems. Numer. Math. 2005, 102, 1-38.

29. Huang, J. Another Version of SOR Iteration and Its Generalization. Master's Thesis, Dalian University of Technology, Dalian, China, June 2013.

30. Li, Y. Study and Realization on Construction, Fairing and Smooth Joining Between Adjacent of the NURBS Surface. Master's Thesis, Xi'an University of Technology, Xi'an, China, March 2010.

31. Surhone, L.; Timpledon, M.; Marseken, S. Non-Uniform Rational B-Spline; Betascript Publishing: Montana, MT, USA, 2010.

32. Zhang, H.; Jiang, D.; Ding, Y. A weight-based optimal faring algorithm for planar cubic NURBS curves. J. Southwest Minzu Univ. 2005, 31, 351-355.

33. Li, X.; Wu, Y.; Leng, H. Research on a New S-curve Acceleration and Deceleration Control Method. Modul. Mach. Tool. Autom. Manuf. Technol. 2007, 50-53.