## Abstract

In spline toolpath interpolation, a crucial point is solving the mapping between
the spline parameter (u) and actual arc length (s) accurately, so that the
toolpath is traveled without undesirable fluctuations or discontinuities in the
feedrate profile. To achieve this, various techniques have been proposed in
literature, including Taylor series interpolation, iterative numerical methods,
and approximating the mapping between u and s with a feed correction polynomial.
This paper presents a new way to parameterize the seventh order feed correction
polynomial, which was introduced by Erkorkmaz and Altintas (2005, “Quintic
Spline Interpolation With Minimal Feed Fluctuation,” ASME J. Manuf. Sci. Eng., **127**(2), pp. 339–349). The proposed technique has a closed-form
solution that can be efficiently implemented in real-time, rather than having to
construct and solve a linear equation system with 14 unknowns for each spline
segment. In this paper, the new solution is derived step by step, and simulation
case studies are presented which demonstrate that the new method accurately
parameterizes the feed correction polynomial in approximately 43% less
computational time, compared to applying the former solution of Erkorkmaz and
Altintas (2005, “Quintic Spline Interpolation With Minimal Feed Fluctuation,”
ASME J. Manuf. Sci. Eng., **127**(2), pp. 339–349). This is because
matrix multiplication operations and a dedicated linear equation solver, which
are cumbersome to implement inside a real-time computer numerical controller
(CNC), are avoided in the new solution.

## Introduction

Spline toolpaths offer many advantages in machining freeform shaped components with complex geometry, such as those produced in aerospace, die and mould, automotive, and biomedical applications. In addition to subtractive machining operations [1], they can also be applied to additive manufacturing [2]. Compared to using linear or circular segments, spline toolpaths enable smoother geometry, better surface finish, and also shorter cycle time due to allowing more aggressive feeds and accelerations to be specified. However, spline toolpaths present a challenge in interpolation, due to the inherent discrepancy between the variable ($u$), which they are parameterized with respect to, and the actual arc length ($s$) that is traveled, as the spline parameter $u$ is incremented.

Figure 1 shows a quintic (fifth degree) spline
toolpath, fit in *x* and *y* axes, formulated in the
form

*k*th segment. As the spline parameter $u$ is incremented, this results in an arc displacement $s$ along the toolpath. The differential relationship between these two variables is

Earlier research in spline toolpath generation has shown that for polynomial type splines (e.g., cubic, quintic, etc.), there exists no known method that can parameterize $u$ to be equivalent to $s$ throughout the whole toolpath. However, polynomial splines are preferred in industrial CNCs due to their simplicity in formulation and implementation. Hence, when the spline parameter $u$ is interpolated in uniform increments, the discrepancy between $u$ and $s$ exhibits itself in the form of two detrimental effects [3]:

- (i)
feedrate fluctuations within each segment;

- (ii)
feedrate discontinuity at segment the connections, which leads to acceleration and jerk spikes.

To alleviate these problems, several approaches have been proposed in literature.

One approach has focused on improved parameterization methods for polynomial splines,
which has led to the nearly arc length parameterized [4], approximately arc length parameterized *C*^{3} [5], and
optimally arc length parameterized [3] quintic
splines. These methods have shown significant reduction in the feed fluctuation and
discontinuity, compared to traditional cubic splines that are currently used in
industry. However, they require extra off-line computation.

Another approach has explored different spline structures, such as
Pythagorean-hodograph curves [6,7], which allow the arc length to be
analytically computed, thereby avoiding the feed inconsistency problem altogether.
However, there is still a strong trend in industry to stick with cubic and
polynomial splines, which are most easily understood, parameterized, and implemented
by practicing engineers. Furthermore, cubic splines provide the necessary degree of
freedom to represent toolpaths with *C*^{2} continuity, which
is sufficient for many of the precision motion control applications in industry.

To improve the feedrate accuracy and continuity during the interpolation of preparameterized toolpaths (such as cubic B-splines), improved interpolation methods have also been proposed. Figure 2 illustrates three common techniques for spline interpolation. The first one is natural interpolation, which is the easiest approach, but it leads to the mentioned feed fluctuation and discontinuity problems. The second one is the approximation of the spline parameter ($ui$), by expanding it into a Taylor series as a function of time, around its value used in the previous interpolation cycle ($ui-1$). This approach, which has been proposed in first order [8–10] and second order [11] forms, unfortunately suffers from accumulating round-off errors. These can be problematic when interpolating long toolpaths, unless adequate numerical measures are incorporated into the practical implementation. The third approach, shown in Fig. 2(c), approximates the nonlinear relationship between $u$ and $s$ with the so-called “feed correction polynomial,” as proposed in Ref. [3]. Typically, seventh degree is used, which allows the six boundary conditions on $u$, $us$ ($=du/ds$), and $uss$ ($=d2u/ds2$) to be imposed at the segment connections; thus avoiding the velocity and acceleration discontinuity problem at these points. The additional two degrees of freedom help establish a reasonably close fit between the feed correction polynomial and the spline parameter ($u$) versus arc length ($s$) data, which is generated earlier on during the numerical integration of the total arc length ($Sk$) for each spline segment. Fitting of the feed correction polynomial is shown in Fig. 3.

In recent research, the feed correction polynomial has also been applied in an adaptive manner to nonuniform rational B-splines (NURBS) interpolation, in order to constrain the mean squared value of the approximation error for the chord parameter [12]. It has also been extended to the ninth degree, in order to interpolate quintic B-splines with jerk continuity [13]. Other polynomial approximations between the spline parameter and arc length have also been proposed in literature, such as the remapping scheme which employs cubic Hermite splines, as reported by Lei et al. [14].

Aside from the above mentioned solutions, iterative numerical techniques have also been explored in Refs. [3,15]. These techniques yield the most accurate results, but typically require at least 2–3 times additional computational time, compared to natural or feed correction type interpolation. Additional measures also need to be incorporated to ensure the stability of iterative algorithms, especially those relying on the Newton–Raphson method. For brevity, in this paper iterative techniques are not illustrated or benchmarked. It was also demonstrated in Ref. [3] that using the feed correction polynomial can provide a good initial guess for a Newton–Raphson based iterative solution.

Overall, feedrate fluctuation avoidance methods, whether they are based on improved toolpath parameterization, Taylor series approximation of $u$ and $s$, polynomial function mapping, or real-time iterative solution, target accurate achievement of the desired arc displacement ($s$), and therefore the correct feedrate $ds/dt$ at a given point along the fitted toolpath. These methods are meant to be used in conjunction with feedrate planning algorithms (or so-called look-ahead functions), like those reported in Refs. [12,14], which optimize the feed profile along the toolpath typically subject to axis velocity, acceleration, and jerk constraints. While feedrate planning is a much broader subject, outside the scope of this technical brief, without accurate interpolation methods, such as the feed correction polynomial approach adopted in this paper, the full advantages of careful jerk and acceleration limiting cannot be appreciated in practice.

In the remainder of this technical note, brief formulation of the feed correction polynomial is reviewed in Sec. 2. Only the essential information, required to follow the proposed new formulation is presented, and details of secondary importance already available in Ref. [3] are skipped. Section 3 presents the derivation of the new method for computing the feed correction polynomial. Section 4 presents implementation results, which confirm the correct solution and effectiveness of the feed correction polynomial, and provide a computational load and time benchmark with respect to the earlier method in Ref. [3].

## Background on Fitting of the Feed Correction Polynomial

Above, $sj$ ($j$ = 1,…, $Mk$) is the numerically integrated arc displacement corresponding to the spline parameter value $uj$. $Mk$ is the number of divisions, which are distributed uniformly along [0,$uMk$]. For each spline parameter value $uj=j\xb7\Delta u$, the corresponding increments in axis positions ($\Delta xj$,$\Delta yj$) are calculated using the toolpath definition in Eq. (1) and the resulting arc displacement is calculated as $sj=sj-1+\Delta xj2+\Delta yj2$, leading up to $sMk=Sk$. $Mk$ is chosen so that the segment with the shortest parameter range (e.g., chord distance) is split into a minimum number of divisions (e.g., $Mk,min$ = 100), and segments with longer chord length are split into proportionally larger numbers of divisions.

While minimizing $J$ in Eq. (6), it is also important to impose the correct boundary conditions for $u$, $us$, and $uss$. This is needed in order to maintain the continuity of the axis velocity and acceleration profiles at the transition points between adjacent toolpath segments $k$ and $k+1$, which are to be interpolated with their respective feed correction polynomials.

where $r(u)=[x(u),y(u)]T$ denotes the position vector. $ru$ and $ruu$ are the derivatives with respect to the spline parameter $u$, which can be
obtained by differentiating, e.g., Eq. (1). Subscripts $k$ and $k+1$,
in the notation $[\u2026]k$, $[\u2026]k+1$,
designate the *k*th and *k* + 1st toolpath segments,
respectively. Thus, it can be validated that the toolpaths are second order
continuous (*C*^{2}) with respect to the spline parameter,
both within each individual segment, and also at the connection points between the
adjacent segments.

Considering Eq. (8), if the commanded feedrate (i.e., tangential velocity $s\xb7$) and tangential acceleration ($s\xb7\xb7$) are continuous throughout the whole toolpath, then the required condition for axis level velocity and acceleration profiles to remain continuous is that the derivatives $us$ and $uss$ also remain continuous.

The coefficients $p0$,…, $p8$ are computed by substituting the expressions for $xu=dx/du$, $yu=dy/du$, $xuu=d2x/du2$, $yuu=d2y/du2$ from Eq. (1) and collecting the terms that multiply different powers of $u$. If the toolpath is cubic, then $P(u)$ is only fourth degree.

Examining Eqs. (9) and (10), it can be seen that due to the conditions in Eq. (7) being satisfied, when a toolpath is interpolated by solving the exact value of the spline parameter ($u$) which yields the desired arc displacement ($s$), and the commanded $s\xb7$ and $s\xb7\xb7$ are continuous, the axis velocity and acceleration profiles in Eq. (8) will always remain continuous, even at the segment crossings. However, such an interpolation calls for a very accurate solution technique, such as the iterative solutions reported in Refs. [3,15].

In the practical implementation in Ref. [3], to avoid numerical conditioning problems, the expressions for $u$ and $s$ were normalized before the computations, and the spline coefficients denormalized afterward. In this paper, the same computer code has been reused for the benchmark study presented in Sec. 4. The normalization steps, for brevity reasons, have not been shown in the preceding review of the earlier method for fitting the feed correction polynomial.

## Proposed New Solution

zeroth, first, and second derivative boundary conditions are already known, as explained in Sec. 2. Hence, solving $usssinit$ and $usssfinal$ would enable the feed correction polynomial to be completely parameterized.

*j*th element of vector $Y$, defined in Eq. (18) as $yj=uj-u\u02dc(\sigma j)$. Thus, $usssinit$ and $usssfinal$ are solved per Eq. (20), which is constructed with the “$\Sigma $” terms calculated from Eq. (19)

Upon solving $usssinit$ and $usssfinal$, the feed correction polynomial coefficients $Akf$,…, $Hkf$ are calculated using Eq. (15).

## Implementation Results

The proceeding implementation results present the following:

- (1)
correctness of the formulation and implementation, by comparing the result with that of the earlier solution in Ref. [3], using a fan shaped toolpath which had also been used in Ref. [3]

- (2)
justification of using a seventh order feed correction polynomial, as opposed to a lower order (such as fifth)

- (3)
arithmetic computational load and computational time comparison with the earlier solution in Ref. [3]

- (4)
an additional toolpath example (based on an airfoil shape), which further demonstrates the effectiveness of the proposed streamlined solution

### Numerical Verification of the Proposed Strategy.

The proposed fitting algorithm for the feed correction polynomial was implemented in matlab and benchmarked against some of the earlier methods that were compared in Ref. [3]. For consistency, the comparison was carried out using the well-known clover leaf profile, originally used by Wang et al. [4,5]. This toolpath, shown in Fig. 4(a), had been used to compare different interpolation algorithms in Ref. [3], and provides a suitable benchmark as it contains a good mixture of low and high curvature regions. Rather than applying feedrate modulation, a constant feed of 100 mm/s was commanded, which makes it very easy to spot and compare feedrate fluctuations consistently throughout the toolpath in terms of percentage (%) values. Implementation results for natural interpolation, first order Taylor interpolation, feed correction based interpolation (using the earlier fitting algorithm in Ref. [3]), and feed correction based interpolation (using the new formulation) are presented in Figs. 4(b)–4(e), respectively. All results were obtained using 1 (ms) as the interpolation period. The improvements obtained with Taylor and feed correction based interpolation are clear. Furthermore, it is seen that the new formulation generates virtually identical results with the earlier method in Ref. [3], which was based on the solution of a 14-unknown linear equation system for each toolpath segment. These were solved using the backslash “\” operator in matlab for solving linear equation systems. Figure 4(e) essentially demonstrates the correctness of the new formulation presented in Sec. 3.

### Justification of Seventh Order Feed Correction.

As indicated in Sec. 1, seventh degree feed correction polynomial was chosen to guarantee acceleration level continuity in the interpolated spline toolpath, while also successfully approximating the mapping between the spline parameter $u$ and actual arc displacement $s$.

The effect of reducing the feed correction polynomial order to a lower degree (i.e., such as fifth order) is demonstrated in Fig. 5. Figure 5(a) shows the feed consistency achieved with the seventh degree case for the fan profile in Fig. 4(a). Figure 5(b) shows the case in which a fifth degree feed correction polynomial is adopted. Of the six degrees of freedom available, four of them were used to match zeroth and first order boundary conditions at the segment connections, and the remaining two to obtain a least squares fit to the generated $u$ versus $s$ data. The degradation from the seventh order case in Fig. 5(a) is clear. Figure 5(c) shows another implementation of fifth order feed correction, this time in which all six degrees of freedom are used only to match zeroth, first, and second derivative boundary conditions at both ends of each spline segment, thus leaving no extra degree of freedom available for curve-fitting to the generated $u$ versus $s$ data. In this case, the degradation in feedrate consistency is even more obvious. These results demonstrate the justification of using a seventh order feed correction polynomial for acceleration-continuous spline interpolation with minimal feedrate fluctuation. If jerk continuity is also desired, a higher order feed correction polynomial, such as ninth, may also be used as proposed in Ref. [13]. The solution methodology for on-the-fly fitting, explained in Sec. 3, can be directly extended to the ninth order feed correction polynomial case as well.

### Computational Load Analysis and Benchmarking.

Arithmetic operation breakdowns comparing the earlier feed correction polynomial fitting method [3] and the proposed new technique are provided in Tables 1 and 2. The analysis in Table 1 assumes that the system of 14 linear equations in Eq. (14) is solved using lower-upper (LU) decomposition with partial pivoting (i.e., Crout's method) [16]. This is one of the mainstream methods available in numerical analysis, which suits the problem and can also be implemented efficiently in real-time. As mentioned earlier, $Mk$ is the number of divisions in arc length integration, during which $u$ versus $s$ data is computed. Using the method in Ref. [3], fitting the feed correction polynomial for one spline segment requires $144Mk+2328$ arithmetic operations. With the new formulation, the corresponding number of operations has been reduced to $84Mk+7$. As $Mk$ is typically is chosen to be larger than 100, the new formulation predicts 41–49% reduction in the number of calculations required. It is important to note that the reasons behind this computational load reduction are: (i) the new streamlined solution does not require the use of a dedicated linear equation solver and (ii) it contains significant simplifications which avoid the tedious matrix multiplication steps, earlier required to construct the $\Phi T\Phi $ and $\Phi Tu$ terms in Eq. (14).

Equation no. | Term | Multiplications | Divisions | Additions | Subtractions |
---|---|---|---|---|---|

(5) and (14) | $\Phi T\Phi $ | 64 × $(Mk+1)$ | 0 | 64 × $Mk$ | 0 |

(12) | $Sk2\u2003\u2026\u2003$$Sk7$ | 6 | 0 | 0 | 0 |

(12) | $L$ or $LT$ | 11 | 0 | 0 | 0 |

(5) and (14) | $\Phi Tu$ | 8 × $(Mk+1)$ | 0 | 8 × $Mk$ | 0 |

(14) | $[\Phi T\Phi LTL0][\theta \Lambda ]=[\Phi Tu\xi ]$ | 1197 | 41 | 0 | 1001 |

Total operations | $72Mk+1286$ | 41 | $72Mk$ | 1001 |

Equation no. | Term | Multiplications | Divisions | Additions | Subtractions |
---|---|---|---|---|---|

(5) and (14) | $\Phi T\Phi $ | 64 × $(Mk+1)$ | 0 | 64 × $Mk$ | 0 |

(12) | $Sk2\u2003\u2026\u2003$$Sk7$ | 6 | 0 | 0 | 0 |

(12) | $L$ or $LT$ | 11 | 0 | 0 | 0 |

(5) and (14) | $\Phi Tu$ | 8 × $(Mk+1)$ | 0 | 8 × $Mk$ | 0 |

(14) | $[\Phi T\Phi LTL0][\theta \Lambda ]=[\Phi Tu\xi ]$ | 1197 | 41 | 0 | 1001 |

Total operations | $72Mk+1286$ | 41 | $72Mk$ | 1001 |

Equation no. | Term | Multiplications | Divisions | Additions | Subtractions |
---|---|---|---|---|---|

(15) | $1/Sk$…$1/Sk7$ | 6 | 1 | 0 | 0 |

(16) | $\sigma j$ | ($Mk-1$) | 0 | 0 | 0 |

(17) | $Sk2$,$Sk3$ | 2 | 0 | 0 | 0 |

(17) | $\sigma j2$… $\sigma j7$ | 6 × ($Mk-1$) | 0 | 0 | 0 |

(17) | $\phi i(\sigma j)$ | 4 × ($Mk-1$) | 0 | 3 × ($Mk-1$) | 1 × ($Mk-1$) |

(17) | $\phi f(\sigma j)$ | 3 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 2 × ($Mk-1$) |

(17) | $u\u02dc(\sigma j)$ | 28 × ($Mk-1$) | 0 | 13 × ($Mk-1$) | 11 × ($Mk-1$) |

(18) | $yj$ | 0 | 0 | 0 | 1 × ($Mk-1$) |

(19) | $\Sigma ii$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma if$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma ff$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma iy$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma fy$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(20) | $usssinit$,$usssfinal$ | 6 | 1 | 0 | 3 |

(15) | $Akf$…$Hkf$ | 42 | 0 | 18 | 12 |

Total operations | $47Mk+9$ | 2 | $22Mk-4$ | $15Mk$ |

Equation no. | Term | Multiplications | Divisions | Additions | Subtractions |
---|---|---|---|---|---|

(15) | $1/Sk$…$1/Sk7$ | 6 | 1 | 0 | 0 |

(16) | $\sigma j$ | ($Mk-1$) | 0 | 0 | 0 |

(17) | $Sk2$,$Sk3$ | 2 | 0 | 0 | 0 |

(17) | $\sigma j2$… $\sigma j7$ | 6 × ($Mk-1$) | 0 | 0 | 0 |

(17) | $\phi i(\sigma j)$ | 4 × ($Mk-1$) | 0 | 3 × ($Mk-1$) | 1 × ($Mk-1$) |

(17) | $\phi f(\sigma j)$ | 3 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 2 × ($Mk-1$) |

(17) | $u\u02dc(\sigma j)$ | 28 × ($Mk-1$) | 0 | 13 × ($Mk-1$) | 11 × ($Mk-1$) |

(18) | $yj$ | 0 | 0 | 0 | 1 × ($Mk-1$) |

(19) | $\Sigma ii$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma if$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma ff$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma iy$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(19) | $\Sigma fy$ | 1 × ($Mk-1$) | 0 | 1 × ($Mk-1$) | 0 |

(20) | $usssinit$,$usssfinal$ | 6 | 1 | 0 | 3 |

(15) | $Akf$…$Hkf$ | 42 | 0 | 18 | 12 |

Total operations | $47Mk+9$ | 2 | $22Mk-4$ | $15Mk$ |

While computing the results reported for the fan profile in Fig. 4, the minimum number of divisions, corresponding to the shortest spline segment, was chosen as $Mk,min=100$. As a result, the average value of the number of integration divisions per segment was $Mk\u2245440$. Based on this number, the total number of arithmetic operations required to implement the earlier feed correction fitting method in Ref. [3], versus the proposed new technique, has been tabulated in bar graph form in Fig. 6(a). In this figure, the abbreviations “M,” “D,” “A,” and “S” designate the required number of multiplication, division, addition, and subtraction operations, respectively. As can be seen, in fitting the 88 feed correction polynomials for the fan profile, the earlier solution requires on average 65,688 operations per segment, whereas the new solution requires 36,967 operations; thus achieving an overall 43.7% reduction in the computational load.

To validate this prediction, computational time benchmarks were carried out for the cases reported in Figs. 4(d) and 4(e). A “for” loop was incorporated into the code, requiring each segment's feed correction polynomial to be solved using both the old and new methods 10,000 times, in order to allow averaging of the measured computational times. The implementation platform was a 2.26 GHz dual core Intel CPU, running matlab r2011a over windows 7 (64-bit) operating system. While implementing the earlier solution in Ref. [3], two different approaches were tested for solving the system of 14 linear equations in Eq. (14). The first approach was using matlab's built-in implementation of LU decomposition with partial pivoting (i.e., Crout's method), by invoking the command “linsolve.” The arithmetic complexity of this method was analyzed in Table 1 and Fig. 6(a). The second approach was to use matlab's backslash \ operator. The linear equation system in Eq. (14) can be expressed as $Ax=b$. According to matlab's documentation [17], when the matrix $A$ is Hermitian but not positive definite (which was numerically verified to be the case for Eq. (14)), matlab applies the block lower triangular diagonal matrix (LDL') factorization method for Hermitian indefinite matrices, available from the MA57 library of routines. While an arithmetic operation breakdown of the LDL’ factorization method is beyond the scope of this paper, computational time benchmarks are still presented in order to allow comparison with the results from Crout's method and the proposed on-the-fly fitting method.

The third case that was benchmarked corresponds to using the new feed correction
polynomial fitting method, proposed in this paper. The average solution times
for all three cases in fitting the feed correction polynomial for a single
segment are shown in Fig. 6(b). It is seen that compared to using the earlier
method in Ref. [3] with Crout's algorithm,
which requires 214 *μ*s per spline segment, the new formulation
allows this duration to be reduced to 123 *μ*s, indicating 42.5%
reduction in the computational time. It is interesting to note that this
reduction is very consistent with the 43.7% reduction in mathematical
operations, as predicted in the analysis in Fig. 6(a). It is also seen that while switching to
the LDL’ algorithm reduces the computational time from 214 to
195 *μ*s (by 8.8%), this reduction is not as significant as
that achieved with the proposed new formulation. This is because both
implementations of the earlier method still require matrix multiplications to
construct the $\Phi T\Phi $ and $\Phi Tu$ terms.
Overall, these benchmarks validate the improved suitability of the new feed
correction polynomial fitting method, over the earlier approach combined with
various linear equation solvers.

### Additional Implementation Example.

Implementation of the proposed feed correction polynomial fitting method on a second toolpath example is shown in Fig. 7. In this case, an airfoil profile is constructed by scaling the point data, which can be found in Ref. [18], with a factor of 6 × and fitting an approximately arc length parameterized $C3$ quintic spline, following the method in Ref. [5]. Two cases are presented: natural interpolation (Fig. 7(b)) and interpolation with a feed correction polynomial, fitted according to the proposed new method (Fig. 7(c)). In natural interpolation, significant feed fluctuation occurs near the highest curvature bend, and reaches 0.5%, whereas the feed correction polynomial contains the fluctuation to within 0.1%; thus, further demonstrating the effectiveness of the proposed scheme.

## Conclusions

This technical note has presented a new way to fit the feed correction polynomial,
which can be used to interpolate along *C*^{2} spline
toolpaths with minimal feedrate fluctuation in an acceleration-continuous manner.
The new approach has a closed-form solution, which is highly suitable for real-time
implementation, as opposed to requiring matrix calculations to construct and solve a
system of 14 linear equations for
each toolpath segment. The equivalence between the old and new solutions has been
demonstrated in terms of achieved feedrate accuracy on a well-known toolpath
benchmark. By conducting a thorough computational load analysis, accompanied by
execution time measurements, it has been demonstrated that the new solution
requires, on average, 43.7% less arithmetic operations and 42.5% less computational
time.

## Acknowledgment

This research was carried out through the support of NSERC granting agency in Canada.