# Hermite Quartic Splines for Smoothing and Sampling a Roughing Curvilinear Spiral Toolpath

Cédric Leroy, Sylvain Lavernhe, Édouard Rivière-Lorphèvre

► To cite this version:

Cédric Leroy, Sylvain Lavernhe, Édouard Rivière-Lorphèvre. Hermite Quartic Splines for Smoothing and Sampling a Roughing Curvilinear Spiral Toolpath. Applied Sciences, 2024, 14 (17), pp.7492. 10.3390/app14177492. hal-04943936

HAL Id: hal-04943936

https://hal.science/hal-04943936v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Hermite Quartic Splines for Smoothing and Sampling a Roughing Curvilinear Spiral Toolpath

Cédric Leroy<fnref n="1" />, Sylvain Lavernhe<fnref n="2" /> and Édouard Rivière-Lorphèvre <fnref n="1" />,*

**Abstract:** From an industrial point of view, the milling of 2.5D cavities is a frequent operation, consuming time and presenting optimization potential, especially through a judicious choice of the tool trajectory. Among the different types of trajectories, some have a general spiral-like aspect and can potentially offer a reduced machining time. They are called curvilinear trajectories and are obtained by interpolation between structure curves, which are the numerical solutions of a partial differential equation. In this case, the machine tool will connect points, and the trajectory will be made up of small segments. While these trajectories exhibit all the necessary qualities on a macroscopic level for rapid tool movement, the tangential discontinuities at a microscopic scale, inherent in the discretization, significantly increase the machining time. This article proposes a method to reparameterize the structure curves of the curvilinear spiral with a set of C<fnref n="2" /> connected Hermit quartic spline patches. This creates a smooth toolpath that can be machined at an average feedrate closer to the programmed one and will, de facto, reduce the machining time. This article shows that the proposed method increases on two representative geometries of cavities and toolpath quality indicators, and reduces the milling time from 10% to 18% as compared to the PDE curvilinear spiral generation method proposed by Bieterman and Sandström. In addition, the proposed method is suitable for any non-convex pocket, with or without island(s).

**Keywords:** pocket milling; roughing; curvilinear toolpath; hermite quartic spline

## 1. Introduction

Pocketing 2.5D cavities, which is common from an industrial point of view (e.g., in aero-nautical structure parts or aerospace electronic crankcase), is a time-consuming operation.

As shown by [1], manufacturing costs are time dependent and typically represent, for dies and molds, 30% of the total price of a piece [2].

In this context, many studies aim to reduce the machining time. The main lines for optimization are enhancing the understanding of cutting phenomena [3–5], improving tool performance [6], selecting the best tools, defining the optimal milling parameters [7–9], and, finally, enhancing the toolpath geometry [10–12].

From this last perspective, the most common approach is to work on the tool trajectory to increase its speed. Multiple studies (e.g., [13,14]) have demonstrated that tangency discontinuities along the trajectory result in increased machining time. These discontinuities necessitate changes in the velocity vector (both in direction and magnitude) along the trajectory, which challenges the machine tool's acceleration and jerk capabilities. As these capacities are limited, the numerical control system of the machine imposes restrictions on these variations, so the tool requires additional time and distance to achieve the programmed feed rate after negotiating corners or halts along the toolpath. Therefore, it is essential to minimize such discontinuities, for which two methods exist:

<figure><img src="image_2.png" /></figure>

### check for updates

**Citation:** Leroy, C.; Lavernhe, S.; Rivière-Lorphèvre, É. Hermite Quartic Splines for Smoothing and Sampling a Roughing Curvilinear Spiral Toolpath. *Appl. Sci.* **2024**, *14*, 7492. https://doi.org/10.3390/app14177492

Academic Editors: Xichun Luo and Abhilash Puthanveettil Madathil

Received: 14 May 2024

Revised: 4 July 2024

Accepted: 22 August 2024

Published: 24 August 2024

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2024 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 (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Machine Design and Production Engineering Lab, UMONS Research Institute for Materials Science and Engineering, University of Mons, Place du Parc 20, 7000 Mons, Belgium; cedric.leroy@umons.ac.be</footnote>

<footnote n="2">² LURPA, ENS Paris-Saclay, Université Paris-Saclay, Avenue des Sciences 4, 91190 Gif-sur-Yvette, France; sylvain.lavernhe@ens-paris-saclay.fr</footnote>

<footnote>* Correspondence: edouard.rivierelorphevre@umons.ac.be</footnote>
* Local adjustments on classical toolpath (Zig-Zag or Contour Parallel) are employed to eliminate or to smooth out the discontinuities [15–17];

* Global strategy to build a spiral-like trajectory (or curvilinear spiral). Macroscopically, this kind of toolpath begins at a center point of the pocket and expands outward without “stop and go”. The toolpath morphs progressively to the shape of the pocket boundary. Natively, curvilinear spirals (or spiral-like toolpath) have a higher degree of continuity [18–21].

Primarily, two approaches are used to build curvilinear spirals:

* Geometrical approaches such as mapping [22,23], morphing [24], or the medial axis method [19,20];

* Differential approaches using an elliptical partial differential equation (PDE) [21] or a level-set propagation model [18]. In those cases, curvilinear trajectories are generated through interpolation between selected structure curves, which are the solutions of the used PDE.

To determine a point on a curvilinear path, an interpolation method is needed. Bieterman et al. [21] use the intersection of a Fresnel vector (winding angle) and two discrete structure curves (Figure 1). The trajectory point is established through a linear angular interpolation between the two intersecting points (Equation (1)):

$$ \rho_S(\theta_i) = \rho_{int}(\theta_i) + (\rho_{ext}(\theta_i) - \rho_{int}(\theta_i)) \cdot \frac{\theta_i}{2\pi}, \quad \theta_i \in [0, 2\pi] \quad (1) $$

with $\rho_S(\theta_i)$, discrete polar equation of a spiral-like turn; $\rho_{ext}(\theta_i)$, discrete polar equation of the external structure curve; $\rho_{int}(\theta_i)$, discrete polar equation of the internal structure curve.

This method is effective when the cavity is sufficiently convex (i.e., without significant protrusions; mathematically, when the polar function of the pocket border is explicit). However, for pockets with more complex geometries, a partitioning task is required as proposed by [25] and illustrated in Figure 2.

Moreover, the structure curves are constructed through interpolations between a series of points calculated using finite element method software. In the context of G1 linear interpolation, the trajectory consists of numerous small segments. Although the trajectory may exhibit characteristics adapted to high feed rates on a macroscopic scale, the tangency discontinuities inherent to this method significantly increase the machining time on a smaller scale. Since ‘polynomial’ interpolators are not yet widely adopted in the industry, efforts are needed to smooth the toolpath and improve its degree of continuity.

<figure><img src="image_4.png" /><figcaption>Figure 1. Interpolation between two structure curves with a Fresnel vector.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Partition of a pocket [25].</figcaption></figure>

The simplest method to achieve this objective is the use of a Finite Impulse Response filter (FIR); [21] suggests a moving average technique. However, the use of such filters may result in a loss of tool position control, potentially leading to a machining in inappropriate zones. Other various studies [22,26,27] achieve three of five-axis trajectory smoothing, via B-splines or Bézier curves. Table 1 shows a summary of the methods proposed in the literature summarizing their main focuses.

<table><caption>**Table 1.** Research timeline.</caption><thead><tr><td>Years</td><td>Authors</td><td>Main Focuses</td><td>Benefits and Limitations</td></tr></thead><tbody><tr><td>&lt;2000</td><td>Held</td><td>Generation of toolpaths on arbitrary geometry Zig-zag or parallel contour</td><td>/</td></tr><tr><td>2003</td><td>Bieterman [21]</td><td>Spiral toolpath PDE approach</td><td>Suitable only for globally convex pocket Not suitable for pocket with island</td></tr><tr><td>2004</td><td>Pateloup [28]</td><td>Corner smoothing</td><td>Local smoothing</td></tr><tr><td>2009</td><td>Held [29]</td><td>Spiral toolpath Geometric approach</td><td>Use of Voronoï diagrams Not suitable for island</td></tr><tr><td>2011</td><td>Xiong [18]</td><td>Spiral toolpath PDE approach Level set propagation</td><td>No smoothing Suitabel for pockets with island but with composite spiral</td></tr><tr><td>2011</td><td>Banerjee [30]</td><td>Spiral toolpath PDE approach</td><td>Only globally convex pocket Suitable for pocket with a centered island G<sup>1</sup> toolpath</td></tr><tr><td>2013</td><td>Xu [23]</td><td>Spiral toolpath Geometric approach Mapping</td><td>Not suitable for pocket with island No smoothing</td></tr><tr><td>2018</td><td>Held [19]</td><td>Spiral toolpath Geometric approach</td><td>Simpler to implement Polygonal spiral (G<sup>0</sup>)</td></tr><tr><td>2023</td><td>Hua [27]</td><td>Curve smoothing B-Spline</td><td>Not integrated in a spiral construction</td></tr></tbody></table>

This article proposes a method to reparameterize the structure curves of the curvilinear spiral with a set of Hermit quartic spline patches C² connected (continuity in tangency and curvature). Analytically, the spiral-like toolpath, formed through linear combination of these curves, maintains the C² continuity, ensuring a smaller machining time. Moreover, the proposed method eliminates the need of Fresnel's vector to establish the spiral parametrization, so the toolpath can be directly created without partitioning the
pocket. It allows a generalization of Bieterman's procedure for any 2.5D pocket geometries, even with non-convex geometry.

## 2. Methodology for Smoothing Toolpaths

Figure 3 summarizes the proposed method and the structure of the following paragraphs.

<figure><img src="image_6.png" /><figcaption>Figure 3. Structure chart of the proposed method.</figcaption></figure>

### 2.1. Determining a Closed Contour

The first step involves establishing a closed contour to define a layer of the 2.5D pocket, which will be milled layer by layer. As the paper is focused on the roughing operation, an offset distance from the border must be accounted for regarding the radial width of the cut of the finishing pass ($D_{offset} = \frac{\phi_{tool}}{2}$).

### 2.2. Solving a Steady State Thermal Problem

Given their smooth nature [31], the solutions of a 2D thermal problem (i.e., isotherms) are chosen for building structure curves: a layer of the pocket is uniformly heated while the border is maintained at 0 [K]. The heat equation is (Equation (2)):

$$ \rho c \frac{\partial T}{\partial \tau} = \nabla(\lambda \nabla T) + q, \quad \partial \Omega = 0 \quad [K] \qquad (2) $$

with $\rho$, the areal density of the layer [kg/m²]; $c$, the specific heat of the layer [J/(kg K)]; $T(x, y)$, the temperature field on the layer [K]; $\tau$, the time [s]; $\lambda$, the thermal conductivity of the material [W/K]; $q$, the areal heat source [W/m²]; $\partial\Omega$, the border condition.

The stationary equation with $q = 1[W/m²]$ and $\lambda = 1[W/K]$ is expressed in Equation (3) and is equivalent to Equation (4) proposed by [21]. Note that the values of $q$ and $\lambda$ do not alter the general shape of the solution temperature field but only impact the maximum temperature:

$$ -\Delta T = \frac{q}{\lambda}, \quad \partial \Omega = 0 \quad [K] \qquad (3) $$

$$ -\Delta\phi(x, y) = 1, \quad \partial\Omega = 0 \qquad (4) $$

### 2.3. Building Smooth Curves from Temperature Field

Equation (3) is solved using a classical finite element method (FEM). Some isotherm curves from the resulting temperature field can be used to construct the spiral-like trajectory. However, due to the discretization inherent in FEM and the necessary interpolations, these curves are microscopically noisy and exhibit tangency discontinuities. To address this issue, a fitting-smoothing method is implemented.
### 2.3.1. Hermite Quartic Spline Interpolation

The Hermite quartic spline interpolation model (HQS) was presented by [32] and is a generalization of the usual Hermite cubic spline interpolation. The general form of the HQS interpolation is (Equation (5)):

$$ \begin{cases} x(u) = ea_0(u)x_0 + ea_1(u)x_1 + eb_0(u) \cos(\theta_0) + eb_1(u) \cos(\theta_1) \\ y(u) = ea_0(u)y_0 + ea_1(u)y_1 + eb_0(u) \sin(\theta_0) + eb_1(u) \sin(\theta_1) \end{cases} \qquad (5) $$

with $P_u = (x(u), y(u))$, the coordinates of an interpolated point; $u \in [0, 1]$, the real inter-polant; $(x_0, y_0)$, the coordinates of the spline first point ($u = 0$); $(x_1, y_1)$, the coordinates of the spline last point ($u = 1$); $\theta_0$ and $\theta_1$, the orientations of the tangent at the beginning and at the end, respectively; and $ea_0, ea_1, eb_0, eb_1$, the blending functions ("e" is for extended).

The degree of blending function equals to 4 and, thereby, 5 constraints are allowed and detailed below.

* The spline passes over two points $P_0$ and $P_1$: +2 constraints;

* The spline respects tangent orientations at the beginning and at the end: +2 constraints;

* An initial curvature $C_0$ is fixed: +1 constraint.

The blending functions are described with four parameters: $k_0, k_1, \alpha, \beta$ (Equation (6), Figure 4). If $\alpha$ and $\beta = 0$, blending functions are those of classical cubic Hermite splines.

<figure><img src="image_7.png" /><figcaption>Figure 4. Representation of blending functions of Hermite quartic splines.</figcaption></figure>

$$ \begin{cases} ea_0(u) = 1 + (\alpha - 3) u^2 + 2(1 - \alpha) u^3 + \alpha u^4 \\ ea_1(u) = (3 - \alpha) u^2 + 2(\alpha - 1) u^3 - \alpha u^4 \\ eb_0(u) = k_0 (t + (\beta - 2) u^2 + (1 - 2\beta) u^3 + \beta u^4) \\ eb_1(u) = k_1 (-\beta + 1) u^2 + (2\beta + 1) u^3 - \beta u^4 \end{cases} \quad u \in [0, 1] \qquad (6) $$

* $k_0$ and $k_1$ do not affect the orientation of the tangent, only the norm of the tangent vector;

* $\alpha$ and $\beta$ are called shape parameters and allow the shape modification of blending functions.
The initial curvature ($C_0$) is a complex function of $k_1, k_2, \alpha$ and $\beta$ (Equation (7)):

$$C_0 = C]_{u=0} = \left. \frac{\frac{dx}{du} \cdot \frac{d^2x}{du^2} - \frac{dy}{du} \cdot \frac{d^2y}{du^2}}{\left( \frac{dx^2}{du} + \frac{dy^2}{du} \right)^{\frac{3}{2}}} \right|_{u=0} = f_0(k_1, k_2, \alpha, \beta) \quad (7)$$

But this expression can be analytically inverted to isolate $\beta$ by the use of the Symbolic Toolbox of MatLab, for example (Equation (8)). The final curvature is obtained in the same way but for $u = 1$ in Equation (7) and is noted as $C_1 = f_1(k_0, k_1, \alpha, \beta)$:

$$\beta = f_0^{-1}(C_0, k_1, k_2, \alpha) \quad (8)$$

### 2.3.2. From Interpolation to Smooth Curve Decomposition

This model of interpolation, between two points, makes it possible to fit and decompose a discrete curve, in this case, an isotherm, into a small number of splines connected in tangency and curvature (named **patches**). Figure 5 explains the procedure followed.

<figure><img src="image_8.png" /><figcaption>Figure 5. General algorithm of isotherm decomposition in minimum $C^2$ connected patches.</figcaption></figure>

$D_{max}$ is the maximum acceptable distance between the curve and patch points; $C_0$ represents the algebraic initial curvature of the first patch. $D_{max}$ and $C_0$ are user-fixed parameters. $D_{opt}$ is the optimized (i.e., minimized) distance between the patch points and HQ spline. **GA** is for the Genetic Algorithm method, and **GM** is for the gradient descent method. $C_1$ is the calculated final curvature of the patch that becomes the initial curvature of the new patch.

By using this procedure, a discrete isotherm may be split into a minimum number of Hermite quartic splines, continuous, connected in tangency and curvature (called **isoHQ**), and respecting a maximum chordal error.

### 2.4. Selecting Structure Curves

To select the structure curves among isotherms and to build the toolpath, isotherms are first decomposed into isoHQ and then sampled.

The analytical expression of patches, along with Equation (9), enables the computation of the curvilinear length of each patch ($s_{patch}$) and, consequently, the total length of an isoHQ. Using normalization (the normalized length is set to one), points can be sampled, equally spaced along an isoHQ:

$$s_{patch} = \int_{0}^{1} \sqrt{\left(\frac{dx(u)}{du}\right)^2 + \left(\frac{dy(u)}{du}\right)^2} du \quad (9)$$

The selection of the structure curves is performed by evaluating the distance ($D_{isoHQ}$) between sampling points of 2 isoHQ (an external one $isoHQ_{ext}$ and an internal one $isoHQ_{int}$)
and linking it with an acceptable distance ($D_{rte}$) equal to the radial tool engagement. If $P_j \in isoHQ_{ext}$ and $Q_i \in isoHQ_{int}$, Equation (10) is computed:

$$D_{isoHQ} = \max_j \{ \min_i [ \sqrt{ (P_j(x) - Q_i(x))^2 + (P_j(y) - Q_i(y))^2 } ] \} \quad (10)$$

* If $D_{isoHQ} < D_{rte}$, the work is restarted with a more internal isoHQ$_{int}$;

* If $D_{isoHQ} > D_{rte}$, the penultimate evaluated isoHQ$_{int}$ is chosen as the structure curve, and the loop is restarted.

Thereby, the distance between two structure curves is under control and always under $D_{rte}$, which guarantees machining without unmilled areas.

## 2.5. Building Curvilinear Spiral

To build a curvilinear spiral, a linear interpolation between points $(xiso_i(j); yiso_i(j))$ and $(xiso_{i+1}(j); yiso_{i+1}(j))$ of the same normalized curvilinear abscissa of two consecutive structure curves is achieved (Figure 6, Equation (11)):

$$ \begin{cases} xs_i(j) = xiso_i(j) + \frac{j}{n}\{xiso_{i+1}(j) - xiso_i(j)\} \\ ys_i(j) = yiso_i(j) + \frac{j}{n}\{yiso_{i+1}(j) - yiso_i(j)\} \end{cases} \quad (11) $$

with $xs_i(j)$ and $ys_i(j)$, the coordinates of the point $j$ from the turn $i$ of the curvilinear spiral, and $n$, the total number of points along the isotherms.

<figure><img src="image_9.png" /><figcaption>Figure 6. Building the curvilinear spiral between 2 isotherms.</figcaption></figure>

## 3. Toolpath Quality Evaluation

As shown in many research works, notably in [13,14], discontinuities in the tangency and curvature of the toolpath geometry challenge the limited acceleration and jerk of each axis of the numerically controlled machine tool. Therefore, the average speed rate is always lower than the programmed one, and the milling time is longer than expected.

Here, a double objective is pursued: increasing the speed along the toolpath (which can reduce milling time) and mastering the position of the tool along the milling path. For the tool speed evaluation, two methods are employed:

* A local approach: VPOp.

* A global approach: Toolpath absolute curvature value.
VPOp (Velocity Profile Optimization—version 2.6) is licensed software developed by X. Beudaert and et al. [33]. Based on kinematic data and the CNC parameters of a given machine tool, it generates optimized feedrate planning that respects the capabilities of the milling center and thus predicts the actual tool speed and slowdowns for the tool. So, the modelized milling time is close to the real one. VPOp is used locally to analyze the velocity profile along two turns of the curvilinear spiral.

To analyze the entire trajectory, histograms of the toolpath absolute curvature value (|C| in mm⁻¹) are used as well as 95% and 99% quantiles of |C| named **Cut95** and **Cut99**. These last indicators can quickly give, unlike VPOp, the geometrical quality of a toolpath and, consequently, the possibility of NC to respect the programmed feed rate.

A comparison will be made between the “no fitted” toolpath (called **RAW**) and the
fitted one with HQS method (called **HQS**). For another comparison point for milling time,
a filtered toolpath is also computed, called **FIR**, and corresponding to a moving average on
three points applied twenty times (this combination gives results close to the HQS method).
Note that the more the filter is applied, the smoother the trajectory but the more the position
control is lost. So, the FIR toolpath is absolutely not suitable to satisfy our goal.

**4. Case Study 1—BIG Pocket**

4.1. *Fitting and Selecting Structure Curves*

The 2.5D pocket given by M. Bieterman and D. Sandström [21] is considered for this work. Because it presents a large variety of tangency and curvature values, this pocket is a good test for the proposed algorithm.

With the conditions exposed in Section 2.4, a temperature field is computed using
a finite element method (12.10³ nodes). The maximum temperature is around 500[K].
The mesh and some isotherms are shown in Figure 7. Our following analysis considers 100
of them; note that the temperature rises from the outside to the inside of the pocket.

<figure><img src="image_10.png" /><figcaption>Figure 7. Mesh and isotherms on a layer of the cavity.</figcaption></figure>

To fit isotherms with HQS patches, the parameter *D*<sub>chord</sub> must be fixed. A too small value increases the number of patches, and the local smoothing is less efficient. On the other hand, a higher value may cause over-machining of the pocket’s border. *D*<sub>chord</sub> = 0.5 [mm] seems, in this case, to be a good compromise but remains a parameter.

At this stage of the work, our fitting method splits each discrete isotherm (composed
with a maximum of about 2500 points) into a maximum of about 25 C² connected patches.
Figure 8 shows the result of a fit along an isothermal curve. Note that the beginning of the
isotherm is chosen in a low-curvature zone. So, the initial curvature of the first patch is
arbitrarily fixed at 0.1 [mm⁻¹] to respect the curvature orientation. As a reminder, the first
patch and the last one of each isoHQ are only C⁰ connected.
<figure><img src="image_11.png" /><figcaption>Figure 8. Local result of the Hermit quartic spline interpolation on an isothermal curve.</figcaption></figure>

Now, a sample of each isotherm as described above must be performed. In total, 360 points are used to represent an isoHQ (Figure 9).

<figure><img src="image_12.png" /><figcaption>Figure 9. Result of the curvilinear sampling.</figcaption></figure>

For this work, the tool diameter is fixed at 10 [mm] with a radial covering of 25%: the acceptable distance ($D_{rte}$) between two isotherms is 7.5 [mm]. Figure 10a shows the 23 automatically selected structure curves.

<figure><img src="image_13.png" /><figcaption>Figure 10. Structure curves and spiral-like tool path.</figcaption></figure>
## 4.2. Building and Evaluation of Curvilinear Spiral

Finally, each revolving of curvilinear spiral is built by linear interpolation between two points of the same index (i.e., the same normalized curvilinear abscissa) of two consecutive structure curves. Figure 10b shows the curvilinear spiral built with our method.

Locally, the feedrate is simulated with VPOp. Two turns of the curvilinear spiral are chosen in three zones of the pocket (internal, intermediate, and external). The modeled machine tool is a *Mikron UCP710*. Its characteristics are reported in Table 2 [34].

<table><caption>**Table 2.** Kinematic characteristics of the Mikron UCP710 machine tool [34].</caption><thead><tr><th>Mikron UCP710</th><th>V<sub>max</sub> [m/min]</th><th>a<sub>max</sub> [m/s<sup>2</sup>]</th><th>J<sub>max</sub> [m/s<sup>3</sup>]</th></tr></thead><tbody><tr><td>X<sub>m</sub></td><td>30</td><td>2.5</td><td>5</td></tr><tr><td>Y<sub>m</sub></td><td>30</td><td>3</td><td>5</td></tr></tbody></table>

For each zone, three fitting methods are tested (RAW = no fit, FIR = moving average, HQS = Hermite quartic spline interpolation), and the toolpath length and the milling time are collected. The average feed rate is calculated. The programmed feed rate (V<sub>m</sub>) is fixed at 10,000 [mm/min]. Data are summarized in Table 3, and Figure 11 shows the velocity profile as a function of the curvilinear abscissa for the three intermediate toolpaths (black = RAW, blue = HQS, red = FIR). Figure 12 links the velocity profile with the positions of the tool along the curvilinear spiral turns. Slowdowns for points A, B, C, E, F, G are due to the “high” curvature of the toolpath. D and H are, in fact, the beginnings of a curvilinear spiral turn, and there, as already explained, the C<sup>2</sup> connectivity is lost.

<table><caption>**Table 3.** Simulated milling time for 3 fitting methods.</caption><thead><tr><th>Toolpath Location</th><th>Fitting Method</th><th>Toolpath Length [mm]</th><th>Milling Time [s]</th><th>Δt/t [%]</th><th>V̄<sub>m</sub> [mm/min]</th></tr></thead><tbody><tr><td rowspan="3">Internal</td><td>RAW</td><td>346</td><td>19</td><td>/</td><td>1092</td></tr><tr><td>FIR</td><td>345</td><td>15.5</td><td>-18</td><td>1335</td></tr><tr><td>HQS</td><td>328</td><td>15.8</td><td>-17</td><td>1245</td></tr><tr><td rowspan="3">Intermediate</td><td>RAW</td><td>650</td><td>22</td><td>/</td><td>1773</td></tr><tr><td>FIR</td><td>645</td><td>17</td><td>-22</td><td>2276</td></tr><tr><td>HQS</td><td>635</td><td>18.5</td><td>-16</td><td>2059</td></tr><tr><td rowspan="3">External</td><td>RAW</td><td>1383</td><td>32.3</td><td>/</td><td>2569</td></tr><tr><td>FIR</td><td>1358</td><td>29</td><td>-10</td><td>2810</td></tr><tr><td>HQS</td><td>1398</td><td>29</td><td>-10</td><td>2892</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Velocity profile as a function of curvilinear abscissa—intermediate turns.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Positions of the feed rate slowdowns along the path (to be linked with Figure 11).</figcaption></figure>

The results in Table 3 demonstrate that the method proposed in the paper, based on the reconstruction of isothermal curves with a small number of Hermite quartic spline patches, reduces the simulated milling time by 10% to 20% in each zone of the pocket, compared to the RAW toolpath. Although the FIR filter is selected for comparison, the HQS method offers an advantage: it maintains control over the distance between the RAW isotherm and the fitted one, a feature not preserved by the FIR method.

The histogram in Figure 13 and the statistical indicators in Table 4 confirm the improvement of the toolpath: small curvatures are more populated, and both Cut95 and Cut99 values are significantly lower for the HQS method. Moreover, this new method eliminates the need for a Fresnel vector to generate the toolpath. Consequently, the limitation of using the Bieterman's method exclusively for 'globally convex pockets' is no longer necessary. Therefore, the HQS method allows the creation of curvilinear toolpaths similar to those by Bieterman and Sandström but applicable to non-convex cavities as demonstrated in the second case study.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** BIG pocket—histogram of $|Curvature|$ for RAW and HQS entire curvilinear spiral.</figcaption></figure>
<table><caption>**Table 4.** BIG pocket—curvature statistical indicators of entire curvilinear spiral.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Curvature Absolute Value (|C|, mm<sup>-1</sup>)</th></tr><tr><th>RAW</th><th>HQS</th></tr></thead><tbody><tr><td>Cut95</td><td>0.24</td><td>0.14</td></tr><tr><td>Cut99</td><td>1.30</td><td>0.97</td></tr></tbody></table>

## 5. Case Study 2—HUG Pocket

The second test case is a pocket with a similar geometry to the pocket proposed by Huang [35]. The major interest of this geometry is its non-convexity: the winding angle (Fresnel's vector, a polar parametrization) approach proposed by Bieterman [21] is not functional, while the method proposed in this paper can handle such geometries. Figure 14 shows that our isotherms reconstruction method works for building a curvilinear spiral for “globally non-convex” pockets. It is built between 17 isotherms and constructed with around 6500 points. Note that the temperature field has two global maxima symmetrically disposed with the pocket axis: two beginning points for the curvilinear spiral are possible.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Entire toolpath.</figcaption></figure>

### *TOOLPATH QUALITY INDICATORS*

Figure 15 shows a histogram of curvature for the entire curvilinear spiral.

As mentioned earlier, the method proposed by Bieterman is unable to generate a toolpath for this pocket. The comparison will be made using only RAW and HQS structure curves (isotherms). As in the previous case, three isotherms are selected (Figure 16): external (iso3), intermediate (iso10), and internal (iso14). The histograms (Figure 17), statistical indicators (Table 5), and simulated milling time (Table 6) demonstrate the geometric improvement in the toolpath quality. In each case, the curvature distribution is more concentrated, and Cut95 and Cut99 values are significantly lower. Additionally, the simulated milling time along each isotherm decreases by approximately 20%.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** HUG pocket—HQS curvilinear entire spiral.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** HUG pocket—structure curves.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** HUG pocket—histogram of curvature absolute value for Iso10.</figcaption></figure>
<table><caption>**Table 5.** HUG pocket—curvature quality statistical indicators.</caption><thead><tr><th colspan="7">Curvature Absolute Value (|C|, mm<sup>-1</sup>)</th></tr><tr><th colspan="7">Entire Spiral</th></tr></thead><tbody><tr><td>Cut95</td><td colspan="6">0.12</td></tr><tr><td>Cut99</td><td colspan="6">0.72</td></tr><tr><th colspan="2">Iso3</th><th colspan="2">Iso10</th><th colspan="2">Iso14</th></tr><tr><th></th><th>RAW</th><th>HQS</th><th>RAW</th><th>HQS</th><th>RAW</th><th>HQS</th></tr><tr><td>Cut95</td><td>0.38</td><td>0.07</td><td>0.45</td><td>0.07</td><td>0.44</td><td>0.10</td></tr><tr><td>Cut99</td><td>1.07</td><td>0.08</td><td>1.05</td><td>0.13</td><td>1.51</td><td>0.12</td></tr></tbody></table>

<table><caption>**Table 6.** Simulated milling time comparison for 3 structure curves.</caption><thead><tr><th>Toolpath Location</th><th>Fitting Method</th><th>Milling Time [s]</th><th>Δt /t [%]</th></tr></thead><tbody><tr><td rowspan="2">Internal (Iso14)</td><td>RAW</td><td>9</td><td>/</td></tr><tr><td>HQS</td><td>7.3</td><td>-19</td></tr><tr><td rowspan="2">Intermediate (Iso10)</td><td>RAW</td><td>15.2</td><td>/</td></tr><tr><td>HQS</td><td>12.2</td><td>-20</td></tr><tr><td rowspan="2">External (Iso3)</td><td>RAW</td><td>22.3</td><td>/</td></tr><tr><td>HQS</td><td>17</td><td>-24</td></tr></tbody></table>

## 6. Conclusions

This article proposes a new method to construct curvilinear spiral toolpaths from structure curves obtained by the FEM solution of PDE inspired from a thermal problem. To reduce the tangency discontinuities due to discretization, the structure curves of the spiral are rebuilt with Hermit quartic spline patches. This method performs the following:

* Models structure curves for any geometry of pockets thanks to the associated thermal model;

* Guarantees a C² continuity of the structure curves and enhances the continuity level of the curvilinear spiral;

* Creates a toolpath having, globally, a low curvature;

* Allows the user to define a maximum chordal error along the structure curves.

Two representative geometries of pockets are used as examples and show reductions of 10 to 18% of machining time as compared to the state of the art. The proposed method also offers a generalization to non-convex pockets as compared to methods from the literature. In the future, the goal will be the use of Hermite quartic curves directly in a numerical controller that can handle spline interpolations.

**Author Contributions:** Conceptualization, C.L.; methodology, C.L.; software, C.L. and S.L.; validation, C.L.; formal analysis, C.L. and É.R.-L.; investigation, C.L.; writing—original draft preparation, C.L.; writing—review and editing, C.L., S.L. and É.R.-L.; visualization, C.L.; supervision, É.R.-L.; project administration, É.R.-L.; funding acquisition, É.R.-L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We wish to thank the LURPA lab from ENS Paris-Cachan for the collaboration proposition and the VPOp licensing agreement, X. Beudaert for the developed software and S. Lavernhe for its parametrization, Fl. Poulaert for the language correction of the paper.
**Conflicts of Interest:** The authors declare no conflicts of interest.

**Glossary**

<table><thead><tr><th>Abbreviation</th><th>Meaning</th></tr></thead><tbody><tr><td>BIG</td><td>Pocket used by Bieterman [21]</td></tr><tr><td>Cut95</td><td>95% quantile of absolute curvature value</td></tr><tr><td>Cut99</td><td>99% quantile of absolute curvature value</td></tr><tr><td>FEM</td><td>Finite Element Method</td></tr><tr><td>FIR</td><td>Finite Impulse Response = a Filtering Method</td></tr><tr><td>GA</td><td>Genetic Algorithm Method</td></tr><tr><td>GM</td><td>Gradient Descent Method</td></tr><tr><td>HQS</td><td>Hermite Quartic Spline interpolation</td></tr><tr><td>HUG</td><td>Pocket geometry used by Huang [35]</td></tr><tr><td>iso</td><td>Isotherm curve</td></tr><tr><td>isoHQ</td><td>Isotherm fitted with HQS method</td></tr><tr><td>PDE</td><td>Partial differential equation</td></tr><tr><td>RAW</td><td>No Fit used</td></tr><tr><td>VPOp</td><td>Velocity Profile Optimization software [33]</td></tr><tr><th>Nomenclature</th><th>Meaning</th></tr><tr><td>C<sub>0</sub></td><td>Initial curvature of a patch (mm<sup>-1</sup>)</td></tr><tr><td>C<sub>1</sub></td><td>Final curvature of a patch (mm<sup>-1</sup>)</td></tr><tr><td>D<sub>isoHQ</sub></td><td>Distance between 2 fitted isotherms (mm)</td></tr><tr><td>D<sub>max</sub></td><td>Acceptable distance between HQS spline and isotherm points (mm)</td></tr><tr><td>D<sub>offset</sub></td><td>Offset distance (mm)</td></tr><tr><td>D<sub>opt</sub></td><td>Optimized distance between patch and isotherm points (mm)</td></tr><tr><td>D<sub>rte</sub></td><td>Radial tool engagement (mm)</td></tr><tr><td>k<sub>0</sub>, k<sub>1</sub>, α, β</td><td>The four parameters of Hermite Quartic Spline interpolation</td></tr><tr><td>u</td><td>The interpolant of Hermite Quartic Splines</td></tr><tr><td>ea<sub>0</sub>, ea<sub>1</sub>, eb<sub>0</sub>, eb<sub>1</sub></td><td>The four blending functions of Hermite Quartic Splines</td></tr><tr><td>s</td><td>Curvilinear length (mm)</td></tr><tr><td>φ<sub>tool</sub></td><td>Diameter of the tool (mm)</td></tr><tr><td>ρ<sub>S</sub>(θ<sub>i</sub>)</td><td>Discrete polar equation of a spiral-like turn</td></tr><tr><td>ρ<sub>ext</sub>(θ<sub>i</sub>)</td><td>Discrete polar equation of an external isotherm</td></tr><tr><td>ρ<sub>int</sub>(θ<sub>i</sub>)</td><td>Discrete polar equation of an internal isotherm</td></tr></tbody></table>

**References**

1. Tlusty, G. *Manufacturing Process and Equipment*, 1st ed.; Pearson: Hoboken NJ, USA, 1999.

2. Altan, T.; Lilly, B.; Yen, Y.C.; Altan, T. *Manufacturing of Dies and Molds*. CIRP Ann.—Manuf. Technol. **2001**, *50*, 404–422. [CrossRef]

3. M'Saoubi, R.; Chandrasekar, S. *Editorial—Special issue "The cutting of metals: Exploring phenomena at the confluence of materials and extreme mechanics"*. Int. J. Mach. Tools Manuf. **2021**, *168*, 103769. [CrossRef]

4. Zhang, D.; Zhang, X.M.; Nie, G.C.; Yang, Z.Y.; Ding, H. Characterization of material strain and thermal softening effects in the cutting process. Int. J. Mach. Tools Manuf. **2021**, *160*, 103672. [CrossRef]

5. Ping, Z.; Yue, X.; Shuangfeng, H.; Ailing, S.; Baoshun, L.; Xiao, Y. Experiment and simulation on the high-speed milling mechanism of aluminum alloy 7050-T7451. Vacuum **2020**, *182*, 109778. [CrossRef]

6. Aurich, J.C.; Kieren-Ehse, S.; Mayer, T.; Bohley, M.; Kirsch, B. An investigation of the influence of the coating on the tool lifetime and surface quality for ultra-small micro end mills with different diameters. CIRP J. Manuf. Sci. Technol. **2022**, *37*, 92–102. [CrossRef]

7. D'Souza, R.M. On setup level tool sequence selection for 2.5-D pocket machining. Robot. Comput.-Integr. Manuf. **2006**, *22*, 256–266. [CrossRef]

8. Aggarwal, S.; Xirouchakis, P. Selection of optimal cutting conditions for pocket milling using genetic algorithm. Int. J. Adv. Manuf. Technol. **2013**, *66*, 1943–1958. [CrossRef]

9. Wu, L.; Li, C.; Tang, Y.; Yi, Q. Multi-objective Tool Sequence Optimization in 2.5D Pocket CNC Milling for Minimizing Energy Consumption and Machining Cost. Procedia CIRP **2017**, *61*, 529–534. [CrossRef]

10. Bouard, M.; Pateloup, V.; Armand, P. Pocketing toolpath computation using an optimization method. Comput.-Aided Des. **2011**, *43*, 1099–1109. [CrossRef]

11. Patel, D.D.; Lalwani, D.I. Quantitative Comparison of Pocket Geometry and Pocket Decomposition to Obtain Improved Spiral Tool Path: A Novel Approach. J. Manuf. Sci. Eng. **2018**, *139*, 031020. [CrossRef]
12. Romero, P.E.; Dorado, R.; Díaz, F.A.; Rubio, E.M. Influence of Pocket Geometry and Tool Path Strategy in Pocket Milling of UNS A96063 Alloy. *Procedia Eng.* **2013**, *63*, 523–531. [CrossRef]

13. Lavernhe, S.; Tournier, C.; Lartigue, C. Kinematical performance prediction in multi-axis machining for process planning optimization. *Int. J. Adv. Manuf. Technol.* **2008**, *37*, 534–544. [CrossRef]

14. Beudaert, X.; Lavernhe, S.; Tournier, C. Feedrate interpolation with axis jerk constraints on 5-axis (NURBS) and {G1} tool path. *Int. J. Mach. Tools Manuf.* **2012**, *57*, 73–82. [CrossRef]

15. Pateloup, V.; Duc, E.; Ray, P. Bspline approximation of circle arc and straight line for pocket machining. *Comput.-Aided Des.* **2010**, *42*, 817–827. [CrossRef]

16. Altintas, Y.; Tulsyan, S. Prediction of part machining cycle times via virtual CNC. *CIRP Ann.* **2015**, *64*, 361–364. [CrossRef]

17. Herzog, R.; Blanc, P. Optimal G2 Hermite interpolation for 3D curves. *Comput.-Aided Des.* **2019**, *117*, 102752. [CrossRef]

18. Xiong, Z.H.; Zhuang, C.G.; Ding, H. Curvilinear tool path generation for pocket machining. *Proc. Inst. Mech. Eng. Part B J. Eng. Manuf.* **2011**, *225*, 483–495. [CrossRef]

19. Held, M.; de Lorenzo, S. On the generation of spiral-like paths within planar shapes. *J. Comput. Des. Eng.* **2018**, *5*, 348–357. [CrossRef]

20. Abrahamsen, M. Spiral tool paths for high-speed machining of 2D pockets with or without islands. *J. Comput. Des. Eng.* **2019**, *6*, 105–117. [CrossRef]

21. Bieterman, M.B.; Sandstrom, D.R. A Curvilinear Tool-Path Method for Pocket Machining. *J. Manuf. Sci. Eng.* **2003**, *125*, 709–715. [CrossRef]

22. Sun, S.; Altintas, Y. A G3 continuous tool path smoothing method for 5-axis CNC machining. *CIRP J. Manuf. Sci. Technol.* **2021**, *32*, 529–549. [CrossRef]

23. Xu, J.; Sun, Y.; Zhang, X. A mapping-based spiral cutting strategy for pocket machining. *Int. J. Adv. Manuf. Technol.* **2013**, *67*, 2489–2500. [CrossRef]

24. Romero-Carrillo, P.; Torres-Jimenez, E.; Dorado, R.; Díaz-Garrido, F. Analytic construction and analysis of spiral pocketing via linear morphing. *Comput.-Aided Des.* **2015**, *69*, 1–10. [CrossRef]

25. Held, M.; Spielberger, C. Improved Spiral High-Speed Machining of Multiply-Connected Pockets. *Comput.-Aided Des. Appl.* **2014**, *11*, 346–357. [CrossRef]

26. Song, D.N.; Ma, J.W.; Zhong, Y.G.; Yao, J.J. Global smoothing of short line segment toolpaths by control-point-assigning-based geometric smoothing and FIR filtering-based motion smoothing. *Mech. Syst. Signal Process.* **2021**, *160*, 107908. [CrossRef]

27. Hua, L.; Huang, N.; Yi, B.; Zhao, Y.; Zhu, L. Global toolpath smoothing for CNC machining based on B-spline approximation with tool tip position adjustment. *Int. J. Adv. Manuf. Technol.* **2023**, *125*, 3651–3665. [CrossRef]

28. Pateloup, V.; Duc, E.; Ray, P. Corner optimization for pocket machining. *Int. J. Mach. Tools Manuf.* **2004**, *44*, 1343–1353. [CrossRef]

29. Held, M.; Spielberger, C. A smooth spiral tool path for high speed machining of 2D pockets. *Comput.-Aided Des.* **2009**, *41*, 539–550. [CrossRef]

30. Banerjee, A.; Feng, H.Y.; Bordatchev, E.V. Process planning for corner machining based on a looping tool path strategy. *Proc. Inst. Mech. Eng. Part B J. Eng. Manuf.* **2011**, *225*, 1578–1590. [CrossRef]

31. Courant, R.; Hilbert, D. *Methods of Mathematical Physics: Partial Differential Equations*; John Wiley & Sons: New York, NY, USA, 2008.

32. Xie, J.; Liu, X. The EH interpolation spline and its approximation. *Abstr. Appl. Anal.* **2014**, *2014*, 745765. [CrossRef]

33. Beudaert, X. Commande Numérique Ouverte: Interpolation Optimisée Pour L'usinage 5 Axes Grande. Ph.D. Thesis, ENS Cachan, Gif-sur-Yvette, France, 2013.

34. Grandguillaume, L.; Lavernhe, S.; Tournier, C. A tool path patching strategy around singular point in 5-axis ball-end milling. *Int. J. Prod. Res.* **2016**, *54*, 7480–7490. [CrossRef]

35. Huang, N.; Lynn, R.; Kurfess, T. Aggressive Spiral Toolpaths for Pocket Machining Based on Medial Axis Transformation. *J. Manuf. Sci. Eng.* **2017**, *139*, 051011. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.