A Method for Generating Toolpaths in the Manufacturing of Orthosis Molds with a Five-Axis Computer Numerical Control Machine

Abstract

This paper proposes a new algorithm for the automatic generation of toolpaths for machining complex geometric positions, such as molds used in orthosis production. The production of individualized orthoses often requires the use of multi-axis machining systems, such as five-axis machines or industrial robots. Typically, complex and expensive CAD/CAM systems are used to generate toolpaths for these machines, requiring the definition of a machining strategy for each surface. While this approach can achieve a reliable and high-quality machining process, it is very time-consuming and makes it challenging to meet the criteria for rapid production of orthopedic aids. Given that their production is a custom-made process using individual shapes as inputs, the toolpath generation process becomes even more demanding. To address these challenges, this paper proposes an algorithm suitable for the automatic generation of toolpaths for such complex positions. The proposed algorithm has been tested and has proven to be robust and applicable.

1. Introduction

Orthoses are aids used to compensate for functional deficiencies in the human locomotor system caused by trauma, pathological processes, deformities, impaired biomechanics, aging, or other factors [1]. The medical–technical field of orthotics and prosthetics (O&P) deals with the production of orthoses [2]. Prosthetics involve products that serve as artificial replacements for missing body parts, while orthoses assist existing body parts in performing specific functions [3].

The construction and production of individual orthoses represent a complex multidisciplinary problem that requires the cooperation of physicians, patients, and appropriate technical experts with the goal of defining the patient’s needs for an appropriate technical solution [4,5]. In addition to understanding the impact of orthoses on the function of each patient’s locomotor system individually, a significant challenge in orthosis production is the material and time resources required, considering that orthoses are still most commonly produced manually [6,7,8,9,10,11,12]. Even after more than 30 years of digital technology application, the need for human labor remains significant because creating toolpaths for individualized positions is difficult to automate, even with modern CAD/CAM systems.

Consequently, the production of such aids represents a significant material cost for the patient and the entire healthcare system [13]. Therefore, the continuous development of CAD/CAM technology in the production of orthoses is necessary to further reduce production costs. The application of this technology has thus far contributed to increased portability; reduced production complexity; greater dedication to working with patients; simplicity in copying and making the same orthoses; and the ease of implementing corrective measures, upgrading the production process, achieving distributed production, reducing storage and production space, and lowering costs while increasing productivity [14,15,16].

Currently, the machine production of orthoses involves the use of CAD/CAM technologies and CNC machine tools or additive manufacturing machines. The advantage of using additive technologies lies in the direct production of orthoses, i.e., eliminating the need for prior mold making. This significantly reduces the need for physical labor and shortens the time and cost of production, which is crucial for the stability of smaller orthotic laboratories [17,18,19]. However, despite advancements in materials and computer solutions, further improvements are needed to fully meet the requirements of orthosis production using additive technologies, such as mechanical properties of the product, production time, and the cost of producing larger volumes.

On the other hand, the application of CAD/CAM technologies in orthotics most often involves 3D surface scanning to determine the geometric model of the mold [20], followed by the production of the mold from polyurethane foam using a CNC machine. The most commonly used machines for making orthosis molds are three-axis or four-axis milling machines, numerical lathes, or their combination with larger working space [21,22]. General-purpose machines are often used, which reduces their implementation cost, but for making larger molds, they need to be modified, representing a significant technical problem and requiring substantial financial investments [23].

Automation and optimization of each step in this process are important for shortening the time required to make a mold or orthosis, so any technical progress in this direction is welcome. Specific technical requirements need to be addressed in the process of making molds/orthoses using CNC machines. The most important ones include the following:

• Fast and simple production of molds/orthoses;
• Implementation of the complex geometry of the workpiece;
• High surface quality with consistent processing precision without vibrations;
• Minimal manual labor;
• Robustness of applied algorithms for generating toolpaths.

To meet these requirements, five-axis machining centers, industrial robotic arms, and special machine designs are increasingly being used [24,25,26,27], which are primarily found in larger laboratories due to their cost.

For mold making, it is necessary to determine adequate cutting toolpaths, defined as a sequence of positions through which the tool passes during material processing to produce the desired shape or workpiece surface finish. Toolpaths are planned and generated using CAM software, which takes into account the geometric model of the mold, types of toolpaths (linear, curved, or circular), strategy (contour, parallel, or z-level milling), and processing parameters (cutting speed, feed rate, and depth of cut). By applying these parameters and adequate mathematical algorithms, toolpaths are optimized according to criteria such as minimum processing time, surface quality, and collision avoidance.

In this context, some research papers have contributed to the development of algorithms for five-axis freeform surface machining. Loney and Ozsoy [28] explore basic methods for generating toolpaths for freeform surfaces using geometric algorithms for interpolation and approximation. These algorithms use splines and polynomials to define smooth toolpaths, which are useful for creating complex orthosis shapes. Elber and Cohen [29] use NURBS (Non-Uniform Rational B-Splines) for precise modeling and toolpath generation for freeform surfaces. Han and Yang [30] employed isophotes for generating toolpaths, enabling uniform processing of complex surfaces. Ding et al. [31] propose an adaptive iso-planar toolpath generation algorithm that adapts to local surface features for optimal processing. Suresh and Yang [32] developed a method ensuring constant waviness height between adjacent tool passes, achieving uniform surface processing.

Currently, CAM modules such as Autodesk Fusion 360 with PowerMill^® 2024 [33], RhinoCAM MILL^® 2024 [34], and similar software are used in practice for making orthosis molds. The most commonly applied methods for achieving specific toolpath planning strategies in commercial CAM programs, focusing on mold surface quality, speed, and simplicity of generating toolpaths for complex surfaces, include the following:

• Contour and profiling (realized using Delaunay triangulation, Voronoi diagrams, Bézier curves, B-spline curves, NURBS curves);
• Three-Dimensional Raster strategy (Uniform rasterization, Digital Differential Analysis (DDA);
• Adaptive clearing (Dijkstra Algorithm [16], A* (Star) Algorithm [35]);
• Voxel-Based Toolpath Generation (voxelization, space search algorithms such as Breadth-first search algorithm).

Preparing toolpaths in these software solutions requires significant time, as it is necessary to define suitable processing strategies for specific features and surface sizes. This approach is common for serial production but is complicated and impractical for making individual products like orthoses and their molds.

Commercial CAM solutions, although effective, often do not offer enough flexibility to adapt existing algorithms for the specific needs of scientific research or applications involving complex surfaces in orthotics. They typically operate with a pre-defined processing strategy designed for broader use but are not always suitable for geometries that require customization. Additionally, while it is desirable for an orthosis mold to be made very precisely, which the aforementioned solutions enable, the tolerances in industrial orthosis production are higher than in metal machining. This allows for the implementation of solutions that can be easily (without demanding interventions) applied for automated toolpath generation of CNC machines for various mold types or orthoses. As a result, various research efforts are focused on developing and analyzing methods to reduce the time needed for computer modeling and toolpath generation. Although these works aim to generally improve existing toolpath generation methods, they are expected to contribute significantly to the orthotic industry.

One current research project for direct toolpath generation in CNC machining centers uses point cloud data. This method involves denoising, down-sampling, and using an octree for processing the 3D model from scanned data [36]. Applying layered slicing and an iterative approach reduces the time needed for surface reconstruction and toolpath planning.

The point cloud-based toolpath generation method is also being developed for robotic processing. Research [28] has shown that toolpath generation can be improved using point cloud projection algorithms. These algorithms allow the automatic generation of adaptive toolpaths, significantly reducing the time needed for manual modeling and adjustment. This method is particularly useful for creating complex anatomical shapes characteristic of orthoses. In this research, a spiral toolpath generation method has been developed for the direct processing of point clouds in three-axis CNC machines. This technique enables efficient surface creation with minimal tool reversals and smooth transitions between layers, which are crucial for the comfort and functionality of orthoses.

Another approach to address the toolpath planning optimization problem for three-axis machining is introduced in [37]. The method consists of two steps. The first step divides the surface into zones with similar local geometric properties to improve the efficiency of toolpath planning algorithms. This is accomplished using an unsupervised clustering algorithm on a mesh defined by isoparametric curves. The second step applies a toolpath planning algorithm to each zone according to its optimal machining direction. This optimal direction is calculated using black-box optimization software. The final goal is to gradually enhance the formulation of the optimization problem to obtain better results.

A method for planning a space-filling curve (SFC)-type toolpath using the Traveling Salesman Problem is proposed in [38]. This method also involves two steps. First, a set of regular cutter contact points is generated on the input surface. A cutting simulation method is developed to evaluate the scallop error and determine the position of the next cutter contact point in the cross-feed direction. Second, the obtained points are input into a Traveling Salesman Problem solver to determine the optimal linking sequences for the cutter contact points.

A novel framework to create a machine learning-based system for choosing the best toolpath planning strategy for CNC machining (finishing) of complex freeform surfaces directly from the CAD model is proposed in [39]. Three toolpath planning strategies are considered: adaptive planar, iso-scallop, and hybrid. First, a novel toolpath analysis module is presented to evaluate the quality of the toolpath considering surface finish, toolpath length, and smoothness. It is then used to analyze and label a large number of CAD models to create a dataset for supervised learning. Finally, a Convolutional Neural Network (CNN) is designed and trained using this dataset to predict the best toolpath planning strategy.

A method for optimizing tool orientation using reinforcement learning has also been proposed for five-axis machining [40]. This approach allows for better tool adaptation to the workpiece surface, reducing processing time and improving manufacturing quality. The algorithms utilize geometric and kinematic information to optimize toolpaths, which is essential for complex orthosis shapes.

By improving CAD/CAM technologies and integrating them into daily production, it is expected not only to enhance efficiency and precision but also to increase the overall quality of orthopedic aids. In the long term, this could lead to better patient outcomes by reducing the time needed for orthosis adjustment and refinement, thereby increasing the comfort and functionality of the final products.

Based on all of the above, the goal of this research was to develop a new algorithm for generating circular-like conformal contour toolpaths for making orthoses or molds for orthoses using a five-axis CNC machine. Its task is to determine the position and orientation of the cutting tool that allows efficient surface milling within the allowable machining error, avoids collisions, and enables the determination of approach vectors for each toolpath point. The proposed solution has several characteristics that significantly contribute to the development of this field, as listed below:

• Automation and efficiency. The proposed algorithm allows for faster automatic toolpath generation, significantly reducing the need for manual intervention and adjustments. This greatly shortens the toolpath preparation time compared to commercial software solutions, which typically require detailed manual adjustments for each model surface.
  This is particularly important, considering that mold production for orthoses often involves complex surfaces requiring precise and customized toolpaths. The automated adaptation of the toolpath, without the need for manual fine-tuning, enables faster and more efficient operations, especially when dealing with the more intricate geometries of orthoses.
• Voxelization flexibility. By using voxelization of the 3D model, the algorithm offers flexibility in choosing voxel size, allowing a balance between precision and processing efficiency. Smaller voxels increase toolpath precision, while larger voxels reduce processing time.
• Toolpath optimization. The algorithm employs multiple methods for toolpath optimization and surface coverage with iso-contours, ensuring efficient material removal and minimizing tool idle time.
• Parallelization and computational efficiency. Voxelization enables parallel computation, which can significantly speed up the toolpath generation process. This is particularly useful on modern hardware platforms that support parallel processing, such as GPUs.
• Robustness and reliability. The algorithm is stable in toolpath calculations, ensuring consistent quality in mold-making for orthoses. This is a significant advantage, as irregularities in the toolpath can lead to inadequate orthosis production, negatively affecting the comfort and functionality of the aid.
• Flexibility. Compared to commercial CAM software, the proposed algorithm is simpler to apply in the O&P industry. Traditional CAM software solutions often requires detailed adjustments for each model surface, while the proposed algorithm features easy customization of its parameters, thereby enabling the adaptation of each toolpath to the specific requirements of the orthosis surface geometry. This is especially important for achieving optimal results in the production of custom orthopedic aids, where greater flexibility in the adaptation of toolpaths is necessary to ensure process optimization in relation to different materials and machining parameters.
• Surface quality. The toolpath generation process allows for achieving high-quality mold surfaces, which are important for making functional and comfortable orthoses.

Finally, it is also important to emphasize that the proposed algorithm ensures the automated generation of equidistant toolpaths with customizable density through an open methodology that allows for modifications and further development. A detailed description of the method, as well as the results of the work in simulation and operational conditions, is presented below.

2. Materials and Methods

The proposed algorithm is performed in five steps (Figure 1):

• Voxelization of the model. The 3D model is divided into small volumetric units (voxels), creating a grid that facilitates further processing.
• Contour extraction using the 3D Region Growing algorithm. An initial point (seed) is selected in the voxel grid, and the algorithm expands the region by adding neighboring voxels that meet similarity criteria.
• Point reduction using the Ramer–Douglas–Peucker algorithm. This algorithm reduces the number of points in the polygon while retaining its basic shape.
• Toolpath optimization using the Traveling Salesman Problem (TSP) algorithm. The points are connected in a sequence to ensure the shortest possible toolpath between all points.
• Toolpath interpolation. The toolpath is generated using B-spline curves to ensure smooth transitions between points.

The main approach chosen to solve the problem is the voxelization of the model made up of a triangle mesh (Figure 2).

Further processing of the obtained results in the form of a set of connected voxels yields contours from which toolpaths are created in subsequent procedures. The realization of contours from the triangle mesh is the subject of numerous studies. Applying algorithms such as Fast Marching or Heat Method in Geodesics [41,42,43,44] achieves high-quality surface coverage with iso-contours. The reason for this choice in our approach is the need to maximize the coverage of the surface to be machined and the fact that the voxelization results will be used in other functions that are necessary for achieving offset surfaces, tool approach directions, collision detection, and other needs.

The proposed algorithm (method) is intended for the machining process to be performed using a five-axis machine or a robotic arm with five or more DOF. The discretization of the model’s surface for the purpose of creating toolpaths has been the subject of several studies [45,46,47]. In this work, triangle meshes stored in STL files were used, which are the standard format for describing the geometry of surfaces captured by a 3D scanner. The basic information for the design of the orthosis is a 3D record obtained by digitization, and this method operates directly with that data. STL files, originally used in additive manufacturing technologies such as rapid prototyping, are a suitable choice of format due to their simplicity in surface description. Unlike the complex descriptions found in the IGES format, STL records define the surface as a triangle mesh, where each triangle is described by the coordinates of its three vertices and a normal vector. This technology provides an important foundation for CAD/CAM applications.

By tiling (tessellating) surfaces and creating groups of triangles, many established methods for clipping, shading, hidden surface removal, and other applications can be utilized. In this work, the approach where the model’s surface is volumetrically discretized (voxelized) was used, forming a voxel grid. In this space, each voxel can connect to 26 other neighbors, a finite number. The choice of raster size is arbitrary; a denser raster implies a greater number of voxels, which increases computation time, particularly in algorithms where, in subsequent steps of the proposed procedure, each point in the series is organized into a connected contour, filtered, and interpolated. Therefore, it is essential to find the optimal voxel dimension to ensure that the entire computation process is not time-intensive while maintaining the quality of the generated toolpaths.

In the proposed algorithm, the wave by which the region grows propagates from a voxel selected according to a specific criterion toward the surrounding connected voxels in a finite number of directions (in the positive and negative X, Y, and Z directions, as well as diagonals). It expands only to the voxels not yet encompassed. The selection of allowed propagation directions, considering the voxel grid structure, can be 6 (orthogonal 3D directions) or 26 (3D orthogonal and diagonal directions). This connectivity allows for uniform spreading in the mentioned directions for each time segment, and the spreading occurs by one voxel per iteration, creating region boundaries one voxel wide. From each boundary, curves with indexed points are formed, serving as the basis for future toolpaths.

Voxelization, as a method whose results create the basis for generating toolpaths using specialized software modules, offers certain advantages, particularly when used for shapes with complex geometry in CNC machining. The resulting structure is organized, making it relatively simple to find neighbors for each voxel. This organization is suitable for further applications in the aforementioned related algorithms, which are particularly important for collision detection and avoidance, as well as determining tool approach directions to the surface being machined.

The proposed algorithm takes this feature into account, as its implementation forms toolpaths with a width of one voxel. This ensures that region boundaries in each iteration are uniformly spaced, which is desirable for toolpath generation. The points are then organized using the Traveling Salesman algorithm, and a spatial B-spline closed parametric curve is passed through them. This method of forming organized toolpaths allows for their reshaping and easier adaptation to different processing strategies based on specific needs.

This process discretizes the space containing the model and enables the creation of a structure where elements are defined as arbitrarily sized cuboids, spatially connected in a way that simplifies propagation or the creation of a tree of the nearest toolpaths between individual segments, or evaluates the nearest distance between two points (voxel centers) on the surface [48,49,50].

2.1. Three-Dimensional Model Voxelization

The process of voxelizing the triangle mesh involves forming a scalar field in which the points represent the centers of the voxels of a predefined size. The voxel grid is a 3D matrix of size N_X × N_Y × N_Z (denoted as voxel_mesh_), where N_X, N_Y, and N_Z are the number of voxel cells along the x, y, and z axes, respectively. The size of each voxel is determined by the grid resolution Δx, Δy, and Δz. Assuming a space with coordinates ranging from (x_min, y_min, z_min) to (x_max, y_max, z_max), the dimensions of the voxel grid will be as follows:

$$\begin{array}{rr}& N_x=⌈{\displaystyle \frac{x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\Delta x}}⌉\\ & N_y=⌈{\displaystyle \frac{y_{\mathrm{m}\mathrm{a}\mathrm{x}}-y_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\Delta y}}⌉\\ & N_z=⌈{\displaystyle \frac{z_{\mathrm{m}\mathrm{a}\mathrm{x}}-z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\Delta z}}⌉\end{array},$$

(1)

where x_min, y_min, and z_min are the coordinates of the bottom left front corner of the grid. Each voxel is filled based on criteria that depend on the application. For example, if a voxel contains a part of the object, it can be assigned a value of 1; otherwise, it is assigned a value of 0.

The algorithm begins by defining the dimensions of the grid and the size of the voxel cells. Then, a 3D matrix voxel_mesh_ is created, in which the values of all voxels are initially set to 0. After that, for each point (x, y, z) in the space, the voxel grid indexes are calculated as follows:

$$\begin{array}{rr}& i=⌊{\displaystyle \frac{x-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\Delta x}}⌋\\ & j=⌊{\displaystyle \frac{y-y_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\Delta y}}⌋\\ & k=⌊{\displaystyle \frac{z-z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\Delta z}}⌋\end{array}.$$

(2)

In the end, the corresponding value is set in the voxel grid for each (i, j, k). The pseudocode (Algorithm 1) for the considered procedure is presented below.

Algorithm 1 Three-Dimensional Model Voxelization—basic and simplified approach

function voxelize(points, voxel_size, grid_min, grid_max):     Nx = ceil((grid_max.x − grid_min.x)/voxel_size.x)     Ny = ceil((grid_max.y − grid_min.y)/voxel_size.y)     Nz = ceil((grid_max.z − grid_min.z)/voxel_size.z)     voxel_mesh_ = zeros(Nx, Ny, Nz)     for each point (x, y, z) in points:       i = floor((x − grid_min.x)/voxel_size.x)       j = floor((y − grid_min.y)/voxel_size.y)       k = floor((z − grid_min.z)/voxel_size.z)       voxel_mesh_[i][j][k] = 1       return voxel_mesh_

This is a basic and simplified formulation of voxelization. However, it can be more complex depending on specific needs, such as interpolation, multiple values per voxel, or different methods of filling and filtering voxel data. In this work, a triangle mesh voxelization approach was used. The dimensions of the voxel grid, in this case, will be as follows:

$$\begin{array}{rr}{N}_x& =⌈{\displaystyle \frac{x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌉\\ N_y& =⌈{\displaystyle \frac{y_{\mathrm{m}\mathrm{a}\mathrm{x}}-y_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌉,\\ N_z& =⌈{\displaystyle \frac{z_{\mathrm{m}\mathrm{a}\mathrm{x}}-z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌉\end{array}$$

(3)

since we chose the basic steps of voxelizing a triangle mesh into a voxel grid with dimensions of ∆x = ∆y = ∆z = 2 mm (2 mm × 2 mm × 2 mm). For each triangle defined by three vertices (x₁, y₁, z₁), (x₂, y₂, z₂), (x₃, y₃, z₃), the boundaries of the voxel cells that are potentially intersected by the triangle can be determined by the following expression:

$$\begin{array}{cc}{x}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{triangle}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(x_1,x_2,x_3\right),& x_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{triangle}}=\mathrm{m}\mathrm{a}\mathrm{x}\left(x_1,x_2,x_3\right)\\ y_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{triangle}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(y_1,y_2,y_3\right),& y_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{triangle}}=\mathrm{m}\mathrm{a}\mathrm{x}\left(y_1,y_2,y_3\right).\\ z_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{triangle}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(z_1,z_2,z_3\right),& z_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{triangle}}=\mathrm{m}\mathrm{a}\mathrm{x}\left(z_1,z_2,z_3\right)\end{array}$$

(4)

After that, the indexes of the voxel grid boundaries can be calculated as follows:

$$\begin{array}{rr}& i_{\mathrm{m}\mathrm{i}\mathrm{n}}=⌊{\displaystyle \frac{x_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{triangle}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌋,i_{\mathrm{m}\mathrm{a}\mathrm{x}}=⌊{\displaystyle \frac{x_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{triangle}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌋\\ & j_{\mathrm{m}\mathrm{i}\mathrm{n}}=⌊{\displaystyle \frac{y_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{triangle}}-y_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌋,j_{\mathrm{m}\mathrm{a}\mathrm{x}}=⌊{\displaystyle \frac{y_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{triangle}}-y_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌋.\\ & k_{\mathrm{m}\mathrm{i}\mathrm{n}}=⌊{\displaystyle \frac{z_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{triangle}}-z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌋,k_{\mathrm{m}\mathrm{a}\mathrm{x}}=⌊{\displaystyle \frac{z_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{triangle}}-z_{\mathrm{m}\mathrm{i}\mathrm{n}}}{2}}⌋\end{array}$$

(5)

All voxel cells within these boundaries need to be checked to determine whether the triangle (from the given triangle mesh) intersects them. If the condition is met, the value of the corresponding voxel cell is set to 1.

Below is an example pseudocode (Algorithm 2) for the described approach. The core part is the function triangle_intersects_voxel, which determines whether a triangle intersects a given voxel cell. This can be implemented using various methods depending on the required precision and execution speed.

Algorithm 2 Three-Dimensional Model Voxelization—a triangle mesh voxelization approach

function voxelize_triangle_mesh(triangles, voxel_size, grid_min, grid_max):     Nx = ceil((grid_max.x − grid_min.x)/voxel_size.x)     Ny = ceil((grid_max.y − grid_min.y)/voxel_size.y)     Nz = ceil((grid_max.z − grid_min.z)/voxel_size.z)     voxel_mesh_ = zeros(Nx, Ny, Nz)     for each triangle (v1, v2, v3) in triangles:       triangle_min = min(v1, v2, v3)       triangle_max = max(v1, v2, v3)       i_min = floor((triangle_min.x − grid_min.x)/voxel_size.x)       i_max = floor((triangle_max.x − grid_min.x)/voxel_size.x)       j_min = floor((triangle_min.y − grid_min.y)/voxel_size.y)       j_max = floor((triangle_max.y − grid_min.y)/voxel_size.y)       k_min = floor((triangle_min.z − grid_min.z)/voxel_size.z)       k_max = floor((triangle_max.z − grid_min.z)/voxel_size.z)       for i in range(i_min, i_max + 1):          for j in range(j_min, j_max + 1):              for k in range(k_min, k_max + 1):                 if triangle_intersects_voxel(triangle, voxel_center(i, j, k, voxel_size)):                    voxel_mesh_[i][j][k] = 1       return voxel_mesh_ function triangle_intersects_voxel(triangle, voxel_center):  # Implementation of an algorithm for checking the intersection of a triangle and a voxel cell  # For example, using SAT (Separating Axis Theorem) or other methods pass function voxel_center(i, j, k, voxel_size):     x = grid_min.x + (i + 0.5) ∗ voxel_size.x     y = grid_min.y + (j + 0.5) ∗ voxel_size.y     z = grid_min.z + (k + 0.5) ∗ voxel_size.z     return (x, y, z)

Assuming there is a 3D voxel grid with dimensions N_X × N_Y × N_Z, each voxel in the grid is represented by a triplet of indexes (i, j, k), where

$$\begin{array}{rr}& {0\le i<N}_x\\ & {0\le j<N}_y\\ & 0\le k<N_z\end{array}.$$

(6)

Common connectivity approaches include the so-called 6-connectivity or 26-connectivity. In basic 6-connectivity, the central voxel being considered can be connected to its neighbors in the −x, +x, −y, +y, −z, and +z directions. In the case of 26-connectivity, diagonal directions at 45 degrees are also included, resulting in a total of 26 neighbors. First, it is checked whether the neighboring voxels are within the bounds of the grid:

$$\begin{array}{rr}& {0\le i+dx<N}_x\\ & {0\le j+dy<N}_y\\ & 0\le k+dz<N_z\end{array}.$$

(7)

For a given voxel (i, j, k), the neighbors include all voxels (i + dx, j + dy, k + dz) where dx, dy, dz take values from the set {−1, 0, 1}, except when dx = dy = dz = 0 (which refers to the voxel itself). An example pseudocode (Algorithm 3) that demonstrates how to find and mark all voxels in contact with the given voxel using 26-connectivity is provided below.

Algorithm 3 Three-Dimensional Model Voxelization—26-connectivity approach

import numpy as np def is_within_bounds(x, y, z, shape):    return 0 <= x < shape[0] and 0 <= y < shape[1] and 0 <= z < shape[2] def get_neighbors(x, y, z, shape)     neighbors = []     for dx in [−1, 0, 1]:        for dy in [−1, 0, 1]:           for dz in [−1, 0, 1]:              if dx == 0 and dy == 0 and dz == 0:                 continue              nx, ny, nz = x + dx, y + dy, z + dz              if is_within_bounds(nx, ny, nz, shape):                 neighbors.append((nx, ny, nz))     return neighbors def voxelize_with_connectivity(voxel_grid):     shape =_voxel_mesh_.shape     _new_mesh = np.zeros(shape, dtype=int)     for x in range(shape[0]):        for y in range(shape[1]):           for z in range(shape[2]):              if_voxel_mesh_[x, y, z] == 1:                _new_mesh [x, y, z] = 1                neighbors = get_neighbors(x, y, z, shape)                for nx, ny, nz in neighbors:                   _new_mesh [nx, ny, nz] = 1     return _new_mesh # Example of a voxel grid (3 × 3 × 3) with one filled voxel in the cente voxel_grid = np.zeros((3, 3, 3), dtype = int) voxel_grid[1, 1, 1] = 1 print("Original voxel mesh:") print(voxel_grid) # Creating a new voxel grid with 26-connectivity _new_mesh = voxelize_with_connectivity(voxel_grid) print("New voxel grid with 26-connectivity:") print(_new_mesh)

The function is_within_bounds checks if the indexes are within the boundaries of the voxel grid and the function get_neighbors returns all neighbors (26-connectivity) for a given voxel (x, y, z) while the function voxelize_with_connectivity iterates through all voxels in the original grid. Moreover, each filled voxel (value 1) must fill the neighboring voxels in the new grid. An example voxel grid is a small 3 × 3 × 3 grid with one filled voxel in the center.

This code creates a new voxel grid where all voxels that are directly or indirectly connected to the initially filled voxel are marked, including all neighbors according to the 26-connectivity rule. In a closed grid, this process defines the shape and volume of the model with greater or lesser precision, depending on the resolution. This feature allows the procedure to be used, for example, for simpler and faster calculations of the moment of inertia for irregularly shaped models, generating toolpaths, finding approach vectors, as well as detecting and predicting collisions, i.e., calculating collisions between machine components and other objects (e.g., cutting tools, machine parts, and other objects in the space). After discretizing the space containing the mesh, only those voxels (cubes) through which the mesh triangles pass are retained. The walls of such a configuration have a thickness of one voxel at any point.

By selecting the voxel size as the fundamental unit, the resolution is defined. Smaller voxels allow for greater precision but increase processing time, as the algorithm must process a larger number of voxels. Conversely, larger voxels reduce the number of computational operations but may compromise the accuracy of the final surface, especially for complex geometries typical of orthoses. To adequately investigate the impact of different voxel sizes, we conducted simulations and compared results with varying voxel resolution settings (1 mm × 1 mm × 1 mm, 2 mm × 2 mm × 2 mm, and 5 mm × 5 mm × 5 mm). The results showed that smaller voxel size (1 mm × 1 mm × 1 mm) provides precise toolpaths but significantly increases processing time, while larger voxels (5 mm × 5 mm × 5 mm) reduce processing time but at the risk of decreased accuracy, which may result in coarser transitions on the final surface.

To achieve a balance between accuracy and speed, a voxel size of 2 mm × 2 mm × 2 mm was selected, as previously mentioned. Preliminary data indicate that with this voxel size, processing time remains within acceptable limits while processing accuracy is maintained at a high level. Medium voxel sizes allow for a sufficiently detailed model to accurately follow complex geometries without excessively prolonging processing time, which could be problematic in industrial applications.

Additionally, the impact of voxel size is also reflected in the surface quality of the workpieces. Smaller voxels enable smoother transitions between toolpaths, reducing the need for additional finishing processes, such as grinding or polishing. Conversely, larger voxels may lead to more visible processing marks, requiring extra time for finishing and negating the time-saving advantage of larger voxels. Regardless of the proposed and selected voxel size in this research, the suggested algorithm provides flexibility in voxel size selection, allowing users to choose the required size based on the specific demands of the project. This flexibility is particularly important in the production of custom orthoses, where different parts of the mold may require varying levels of precision.

2.2. Contour Extraction

To achieve a specific contour passage, this work employs a 3D Region Growing algorithm. This algorithm is used for shape segmentation in models composed of defined volumetric units, with the aim of extracting specific connected structures within such models. The algorithm begins by propagating from a given seed point and expanding the region by adding neighboring voxels that meet certain similarity criteria. In the context of 26-connectivity, each iteration can include only one neighboring voxel per propagation step (Figure 3) [51,52,53,54,55,56,57,58,59].

The process of setting the seed point is arbitrary and is performed manually. Within the simulation software, the model is positioned in the machine’s workspace, and then the seed point is set on the positioned model. The criteria for selecting this point vary depending on the shape and size of the model. In most cases, the point chosen is the one that is Euclidean-farthest from the center of the fixture assembly. To increase system flexibility and enable better tool access to all model surfaces, the model is fixed in the machine’s workspace using only a single fixture assembly, providing more space for free tool movement.

The selection of the current criteria is influenced by the fact that mold models can have diameters exceeding 800 mm and are made from brittle polymer materials. Therefore, the starting point for machining is chosen in such a way that the process begins from the furthest end and moves toward the fixture, where the model is oriented with its wider section. This approach ensures that a sufficiently rigid base is maintained throughout the machining process, preventing breakage of the workpiece due to the forces and moments exerted on it.

The algorithm consists of the following steps [60,61,62,63,64]:

• Initial point: The algorithm starts from the selected seed voxel V_i (i = 0) with coordinates (x₀, y₀, z₀) and intensity I₀. This voxel is added to the candidate set, i.e., the set of voxels that make up the region (R).

  R = {V₀}.

  (8)
• Similarity criterion: The similarity criterion S(V_i, V_j) is defined to determine whether the neighboring voxel V_j should be added to the region. In the proposed algorithm, the intensity is defined as the spatial position of a voxel. Therefore, the similarity criterion is defined as the Euclidean distance between the centers of two voxels—voxel V_i, which is the currently analyzed voxel, and voxel V_j, located in its vicinity:

  S(V_i, V_j) = ∣𝐼(V_i) − 𝐼(V_j)∣.

  (9)
  In other words, the criterion is a distance of 1 voxel from the last voxel in the previous wave, i.e., a distance of 1 voxel from the previously found voxel and a distance of 1 voxel from the voxel in the previous wave.
• Expansion: If, for each voxel V_i in region R and for each of its neighbors V_j, the following condition is met:

  S(V_i, V_j) ≤ T,

  (10)

  where T is the similarity threshold, then voxel V_j is added to region R (V_j ∈ R). In this case, the threshold acts as a tolerance within which two voxel centers can be considered part of the same surface. Since we use the 26-connectivity method, the proposed algorithm checks all 26 neighbors of voxel V_i in 3D space. The similarity criterion is applied to all these neighboring voxels to determine which of them meets the similarity threshold, i.e., geometric proximity to the current voxel. This ensures that the algorithm can track complex shapes and curvatures since each voxel has the potential to connect with neighboring voxels in all directions. For example, if the geometry is highly curved, diagonal voxels can be crucial for maintaining region continuity because direct neighbors (only along the X-, Y-, or Z-axis) may not meet the geometric similarity criterion. This allows the algorithm to accurately expand the region and generates smooth and consistent toolpaths that follow complex surface shapes, which is essential for machining in industrial applications.
• Iteration: The steps are successively repeated until all relevant voxels are included.

The pseudocode of the described 3D Region Growing algorithm is presented below (Algorithm 4).

Algorithm 4 Three-Dimensional Region Growing Algorithm

def region_growing(volume, seed, threshold):     x0, y0, z0 = seed     region = set()     candidates = set([seed])     while candidates:       x, y, z = candidates.pop()       if (x, y, z) in region:             continue       region.add((x, y, z))       for dx in [−1, 0, 1]:          for dy in [−1, 0, 1]:             for dz in [−1, 0, 1]:                if dx == 0 and dy == 0 and dz == 0:                   continue                nx, ny, nz = x + dx, y + dy, z + dz                if (nx, ny, nz) in region:                    continue                if abs(volume[nx, ny, nz] - volume[x, y, z]) <= threshold:                    candidates.add((nx, ny, nz))     return region

If a surface made up of connected voxels were flat, the wave from the starting point would spread concentrically in the form of a discretized circle. However, when the wave propagates across the surface of a structure with a complex shape, the displacement of its spread in a single iteration will be only one voxel. This means that in each phase, the wave’s propagation creates a new edge with a width of one voxel. These edges are connected, and a curve can be drawn through them. This curve is one voxel away from the previously generated curve, allowing for the formation of a set of such curves. This set of curves forms the basis for generating the toolpath. The result of the described procedure is shown in Figure 4, Figure 5 and Figure 6. The colors of the individual contours in Figure 4, Figure 5 and Figure 6 were chosen based on a random sample model among 16 predefined RGB schemes. The reason for selecting different colors for each contour is to visually distinguish each individual curve.

Each curve thus forms a separate toolpath. In the proposed algorithm, due to local geometric features, there is a possibility that wave propagation may lead to merging, which can result in the formation of loops and, subsequently, two or more new curves. In such cases, the algorithm propagates within the area defined by the contours of each loop until all voxels are selected. If, during this process, some curves become “trapped,” the algorithm moves to the next contour, and the process is repeated for each newly generated contour.

At the end of this step, a sequence of unorganized points is obtained (i.e., the order of voxel centers where the XYZ coordinates have integer values). The resulting points may exhibit local deviations in the form of branches that need to be smoothed out to form a uniform sequence without abrupt transitions. Given the structure of the voxel grid, unnecessary points are filtered out in subsequent steps, and the resulting sequence of organized points is further smoothed. In our example, the Ramer–Douglas–Peucker algorithm was used for this purpose [65,66].

2.3. Point Reduction

The Ramer–Douglas–Peucker algorithm reduces the number of points representing a polygon while preserving the basic shape of the line connecting its corner points. The algorithm removes points that are closer to the line than a specified distance threshold. It follows these steps:

• Inputs: A sequence of points P = [P₁, P₂, …, P_n], which represent a polygon, and a distance threshold ε represent input variables to the algorithm.
  Let P₁ = (x₁, y₁, z₁) and P_n = (x_n, y_n, z_n) be the endpoints of the line segment. The distance of the point P_i = (x_i, y_i, z_i) from the line P₁P_n, which connects the first (P₁) and last (P_n) point in space, is given by the following formula:

  $$\begin{array}{c}{d}_i={\displaystyle \frac{\sqrt{A^2+{\mathrm{B}}^2+{\mathrm{C}}^2}}{\sqrt{{\left(x_n-x_1\right)}^2+{\left(y_n-y_1\right)}^2+{\left(z_n-z_1\right)}^2}}},\\ A=\left(y_{\mathrm{i}}-y_1\right)\left(z_n-z_1\right)-\left(z_{\mathrm{i}}-z_1\right)\left(y_n-y_1\right),\\ B=\left(z_{\mathrm{i}}-z_1\right)\left(x_n-x_1\right)-\left(x_i-x_1\right)\left(z_n-z_1\right),\\ C=\left(x_i-x_1\right)\left(y_n-y_1\right)-\left(y_i-y_1\right)\left(x_n-x_1\right).\end{array}$$

  (11)
  The distance of each point from the line P₁P_n needs to be calculated in order to find the point P_k with the maximum distance dmax.
• Iteration (Recursive Decomposition):
  If dmax > ε, retain P_k and recursively apply the algorithm to the two segments: [P₁, P_k] and [P_k, P_n].
  If dmax ≤ ε, remove all points between P₁ and P_n.
  In other words, if dmax > ε, the algorithm is recursively applied for the segment before and after that point. Otherwise, only the endpoints are retained.

By combining the obtained results, a simplified sequence of points is achieved. This algorithm efficiently reduces the number of points while preserving the basic shape of the lines. Its pseudocode is presented below (Algorithm 5).

Algorithm 5 Ramer–Douglas–Peucker Algorithm

function RDP(points, ε):     d_max = 0     index = 0     n = length(points)     # Find the point with the greatest distance     for i from 2 to n − 1:       d = perpendicular_distance(points[i], points[1], points[n])       if d > d_max:           index = i           d_max = d     # If the maximum distance is greater than ε, recursively simplify     if d_max > ε:       # Recursively simplify       rec_results1 = RDP(points[1:index + 1], ε)       rec_results2 = RDP(points[index:n], ε)       # combine results       result = rec_results1[1:end − 1] + rec_results2     else:       result = [points[1], points[n]]     return result function perpendicular_distance(point, line_start, line_end):     # Calculate the distance of a point from a line     numerator = sqrt(((point.y − line_start.y) * (line_end.z − line_start.z) − (point.z −        line_start.z) * (line_end.y − line_start.y))^2 + ((point.z − line_start.z) * (line_end.x −        line_start.x) − (point.x − line_start.x) * (line_start.z − line_start.z))^2 + ((point.x −        line_start.x) * (line_end.y − line_start.y) − (point.y − line_start.y) * ( line_end.x −        line_start.x))^2)     denominator = sqrt((line_end.y − line_start.y)^2 + (line_end.x − line_start.x)^2 +        (line_end.z − line_start.z)^2)     return numerator/denominator

2.4. Toolpath Optimization

After the aforementioned step, the procedure continues by connecting the points of each sequence into an organized (indexed) continuous toolpath. This toolpath ends at the starting point, with each individual point connected solely and exclusively to two other points (Figure 7). To achieve this, the Traveling Salesman Problem (TSP) algorithm was used. Although more efficient methods exist for determining the shortest toolpath among given points, TSP was selected due to its quick implementation for a relatively small number of points, ease of tracking execution, and robustness despite its slower execution.

The TSP is a classic problem in graph theory and optimization, aiming to find the shortest possible route that visits each point exactly once and returns to the starting location. After connecting the points into an organized contour, the nearest point on the next contour is found, and the indices of the points are shifted so that this point becomes the starting point.

• Inputs: A set of “cities”, or in our case, points P = {p₁, p₂, …, p_n} along with the distances between each pair of points d(p_i, p_j) for all i, j where 1 ≤ i, j ≤ n.
• Outputs: A permutation P = (p₁, p₂, …, p_n) of points that minimizes the total traveled distance:

$$\mathrm{m}\mathrm{i}\mathrm{n}\sum _{i=1}^{n-1}d\left(p_i,p_{i+1}\right)+d\left(p_n,p_1\right).$$

(12)

In this algorithm, two common approaches to solving the problem are distinguished: Brute Force and Dynamic Programming. The Brute Force method is the simplest yet least efficient way to solve the TSP, as it tests all possible permutations of cities and calculates the total distances for each permutation. In contrast, a more widely used approach for finding the shortest route among points involves a classic Dynamic Programming solution known as the Held–Karp algorithm. Although this method also has exponential complexity, it allows for more efficient solving of the TSP compared to the Brute Force method, making it significantly faster for a larger number of points [67,68,69,70,71,72,73]. The procedure involves the following steps:

• Define the state: Let dp[S][i] be the minimum cost to visit all cities (points) in subset S ending at city i_i, and i is a city in S (S ⊆ {p₁, p₂, …, p_n}). Variable dp is called the DP table.
• Initialization: For a subset S that contains only one city, you set the value of dp[S][i] to infinity for all i except for the starting city (in which case dp = 0). This is because, initially, the cost to reach any city (except the starting city) from an empty set is not defined (infinity).
• Filling of the DP table: The table is filled in for all subsets using a recursive formula

  $$\mathit{dp}\left[S\right]\left[i\right]=\underset{i\in S,i\ne j}{min}(dp[S\left\{j\right\}\left]\right[i]+d(i,j\left)\right).$$

  (13)
• Solution: At the end, the algorithm finds the minimum distance and reconstructs the optimal route:

  $$\underset{i\ne 1}{min}\left(\mathit{dp}\left[\left\{p_1,p_2,\dots ,p_n\right\}\right]\left[i\right]+d\left(i,1\right)\right).$$

  (14)

The pseudocode for the Held–Karp algorithm is given below (Algorithm 6).

Algorithm 6 Held–Karp Algorithm

import itertools def held_karp(d):     n = len(d)     C = {}     # Initialization     for k in range(1, n):        C[(1 << k, k)] = (d[0][k], 0)     # Filling the dp table     for subset_size in range(2, n):        for subset in itertools.combinations(range(1, n), subset_size):           bits = 0           for bit in subset:              bits |= 1 << bit           for k in subset:              prev_bits = bits & ~(1 << k)              res = []              for m in subset:                 if m == k:                    continue                 res.append((C[(prev_bits, m)][0] + d[m][k], m))              C[(bits, k)] = min(res)     # Finding the minimal route     bits = (2**n − 1) − 1     res = []     for k in range(1, n):       res.append((C[(bits, k)][0] + d[k][0], k))     opt, parent = min(res)     # Reconstructing the route     path = []     for i in range(n − 1):        path.append(parent)        bits = bits & ~(1 << parent)        _, parent = C[(bits, parent)]     path.append(0)     return opt, list(reversed(path))

This step ensures that each sequence (obtained by applying the 3D Region Growing algorithm) is organized in such a way that the total toolpath connecting the points is minimized.

2.5. Toolpath Interpolation

After obtaining a set of unorganized points for each time step (dt) using the 3D Region Growing algorithm, indexing them with the Traveling Salesman algorithm, and filtering local deviations using the Ramer–Douglas–Peucker algorithm, the conditions for implementing the toolpath interpolation procedure are established. This procedure results in smooth toolpaths with uniformly defined segment lengths and allows for the creation of curves with arbitrary resolution and local curvature. The starting positions on each curve are determined to enable the tool to transition smoothly to the next curve without abrupt changes. In this research, toolpath interpolation is performed using a B-spline curve.

In general, the interpolation process is defined by parameters that determine the smoothness of the curve and how it adapts to control points. Given the density of these points, this process is simplified because the B-spline curve is not approximated (which is an optimization process that is computationally expensive and where the position of control points on the curve is not definitively determined). Instead, the centers of individual voxels are directly used as control points for the B-spline, determining the basic shape of the curve. For this procedure to be possible, it is essential to ensure the organization of the points (centers of the voxels) into a continuous sequence without local loops or branches, which must be properly indexed. This allows for the successive joining of a series of points obtained by wave propagation over the voxelized surface, where the start and end of each series are in contact with the same positions of the previous or subsequent series of points.

The selection of the B-spline parametric curve is based on its ability to pass through given points with minimal deviation. Although the curve does not need to pass directly through these points, it describes them well due to the dense distribution of voxels. Knots form a sequence of values (knot vector) that defines how control points are connected. The number of knots depends on the number of control points and the degree of the B-spline curve. B-spline curve parameters also include basis functions, which are defined by the knots. Each basis function contributes to the final curve according to the corresponding control point. The final B-spline curve is obtained by summing all the control points, each multiplied by its corresponding basis function, ensuring the curve is smooth and continuous while adjusting to changes in control points and knots.

In this research, a fourth-order B-spline curve was chosen. This closed-type parametric curve has the characteristic of ensuring C2 continuity, meaning that the first and second derivatives of the curve are continuous over the entire domain. Additionally, one of the leading advantages of B-spline interpolation is its ability to locally adjust the curve without altering its global shape, i.e., by changing the order of the curve and the positions of the control points. All of this is necessary to achieve smooth tool movement during the machining of complex orthosis geometries, which is essential for achieving the required high precision in processing.

The differentiability of the curve enables the generation of local tangent, normal, and binormal vectors, known as the Frenet frame. This frame allows for the local definition of transformations, tool tilt direction, or axis of rotation, which is necessary for the implementation of the proposed algorithm. Using these vectors, it is possible to control the tool’s approach direction and assess the curvature of the surface in specific segments. It is important to note that from the generated toolpaths, compensated toolpaths are formed in subsequent steps. In other words, for each point and its corresponding direction on the final toolpath, whether the tool’s tip can touch the surface is evaluated. The desired direction is the normal vector, initially defined, but due to the surface geometry and potential collisions, this may not always be feasible. The process for determining the tool approach direction is based on identifying the point where the ball end of the tool touches the surface of the model, with the optimal tool center position being iteratively calculated.

The obtained toolpaths do not define the final tool position on the surface, which will be addressed in a subsequent procedure. The toolpaths are intended for five-axis machining, so alongside the realized toolpath points and the calculated tangent, normal, and binormal vectors, alternative approaches to the tool’s interaction with the surface will be considered to avoid potential collisions between the tool, the machine’s moving structures, and the workpiece (Figure 8). For testing purposes, the approach of the tool to the surface using the normal vector of the point will be considered, which is expected to be the most desirable option.

3. Results and Discussion

The proposed algorithm was tested under both simulated and experimental conditions using mold models of the ankle–foot orthosis, orthopedic insoles, and spinal orthosis.

3.1. Simulations

The simulation of the obtained CAD/CAM models of the orthosis mold was conducted using custom-made software (CogniVision^®, ver. 4.5). This program allows for positioning the model with generated toolpaths within the machine’s workspace and inside a previously imported workpiece model (Figure 9 and Figure 10). Upon positioning the model, the program automatically generates a conical shape at the base of the model, which radially expands towards the workpiece holder. This conical form not only reinforces the model but also prevents collisions between the tool and the workpiece holder. Once the necessary parameters are defined, the machining simulation begins.

The simulation procedure includes positioning the workpiece model within the machine’s workspace. The positioning corresponds to its actual location, facilitated by a special feature of the software. Transparency is applied to the workpiece model, allowing for visual monitoring of the positioning and rotation of the mold model within its volume. This step ensures that the mold model is entirely within the dimensions of the oriented workpiece model. Upon completion of the positioning, a conical model is automatically formed in the workpiece holder section of the model.

The toolpath coordinates are linked to the mold model. During manipulation and positioning in the virtual workspace, these points follow the base model and are aligned accordingly within the machine’s workspace. When the main module of the program is activated, the machining simulation is displayed. With the known position and local coordinate system of each point, along with the defined parameters and machine kinematics, the machining process is visualized. This procedure allows for a preliminary assessment of the correctness of the machining process through visual observation. If no errors are detected during the simulation, the procedure is then approved for execution on the actual machine tool.

In the introductory chapter of the paper, one of the key advantages of the algorithm mentioned is its speed in generating toolpaths. To assess this advantage, we conducted an analysis of the time duration for individual phases of the algorithm using a simulation environment with an Intel i7 processor (Intel Corp., Santa Clara, CA, USA), 16 GB of RAM (Kingston Techology Corp., Fountain Valley, CA, USA), and an Nvidia RTX 3080 (Nvidia Corp., Santa Clara, CA, USA) graphics card. Each processing phase was tested on the geometry of a typical orthosis, i.e., the following data were obtained based on the processing of an orthopedic insole model consisting of 281,330 voxels. For this model, the durations of individual phases were as follows:

• Three-Dimensional Model Voxelization—3.45 s;
• Contour Extraction—4.70 s;
• Point Reduction—1.25 s;
• Toolpath Optimization—2.20 s;
• Toolpath Interpolation—3.50 s.

Since these results were achieved without algorithm optimizations or multithreading, the proposed algorithm has the potential for further reduction in toolpath generation time.

3.2. Experimental Work

3.2.1. Testbed Description

After evaluating the functionality of the generated toolpaths in the simulation software, the results of the proposed method were tested on multiple samples using a five-axis machine. A specially designed machine with a parallel kinematic architecture, which enables high speed and precision in the processing of complex surfaces, was developed specifically for the production of molds for orthoses and prostheses and utilized in this research [74].

The machine consists of two main modules: a tripod module based on parallel kinematics and an assembly with rotational axes A and C (Figure 11).

The tripod module includes a base platform shaped as an equilateral triangle positioned at the top of the machine, from which three actuators (linear axes) extend downward and toward the machine’s center. These axes form an inverted pyramidal shape, converging at a support point along the machine’s central axis. The actuator linear axes are set at a 120° angle relative to each other and at a 35° angle relative to the horizontal plane. This geometry affects the kinematic equations and the relationships among the coordinates.

Each linear axis features a slider connected to the moving platform via two arms with spherical joints, and the arms are also connected through spherical joints to the corners of the moving platform. Although the connection between each linear axis and the moving platform is achieved through two arms, functionally, it acts as a connection through a single point. This construction allows the moving platform to move while maintaining a constant distance (parallelism) between the arms, regardless of the position of the moving platform, whose translational motion is caused by changes in the slider positions on the linear guides to which it is connected through fixed-length arms. This design prevents the rotation and tilting of the moving platform, ensuring that its motion is always translational in the X, Y, and Z directions. The plane in which it moves remains parallel to the plane of the base platform throughout the entire operation.

At the base of the machine, there is an assembly with two rotational axes, A and C. These axes are perpendicular to each other, with the plane of rotation serving as the reference surface during the machine’s operation. Axis A is collinear with the machine’s X-axis, which offers certain advantages when solving direct and inverse kinematics. The X-axis of the machine is a line passing through the rotational A-axis when its angle of rotation is 0. The Y-axis of the machine passes through the rotational C-axis when its angle of rotation is set to 0 degrees. The Y-axis is perpendicular to the machine’s X-axis and lies in the same plane. The C-axis features a clamping fixture that enables the securing of a workpiece or other object and its rotation along a defined path. The range of rotation for axis A is from −45° to +45°, while axis C can rotate from −180° to +180°.

The machine’s Z-axis is perpendicular to the plane of the base platform, as well as to the plane in which axes A and C lie when both are in the neutral position. In other words, the Z-axis passes perpendicularly through the intersection of the X and Y axes, which have also been chosen as the origin of the machine’s coordinate system. The tool, mounted on the spindle motor, is centrally located on the moving platform and remains parallel to the machine’s Z-axis at all times. Due to the application of parallel kinematic principles, the movements of the moving platform are achieved by changing the lengths of the actuators.

The linear axes of the machine are controlled via a computer controller that receives instructions in the standard G-CODE format. It is important to note that changes in the lengths of the linear guides do not directly reflect coordinates in the machine’s workspace but instead represent variables in the system. These variables, together with constant parameters, allow for the calculation of the positions of individual points within the workspace. The machine’s kinematics involve a complex nonlinear relationship between the lengths of the actuators (absolute positions of the sliders on the linear axis) and the tool tip position. Direct kinematics allows for determining the tool position from known actuator lengths, while inverse kinematics solves for the required actuator lengths for a given tool position. For kinematic analysis, parameters that describe the platform’s position and movement need to be defined. The main constant parameters that define the geometry and relationships among the axes, enabling the solution of the machine’s inverse and direct kinematics, include the following:

• The spatial positions of the points where the linear axes connect to the perimeter of the base platform in the machine’s coordinate system.
• The direction vectors of the linear axes.
• The current positions of the sliders on the linear axes.
• The constant lengths of the three arms connecting the sliders and the moving platform.
• The positions of the arm connections to the moving platform in the coordinate system of the moving platform, which are centered at the middle point and have their primary axes parallel to the machine’s coordinate system. Since the moving platform does not rotate or tilt during operation, these points’ positions are determined by the position of its center point, onto which the tool tip is orthogonally projected. The tool tip is offset from the platform by the tool’s length and an appropriate offset corresponding to the clamping length and the spindle motor shaft segment length, measured profile-wise from the moving platform plane to the tool clamp connection.
• The reduction in tool length by the radius of its tip.
• The radius of the tool tip—in our case, a sphere since a ball nose end tool is used.

3.2.2. Machining Parameters

A cylindrical ball nose end mill with a diameter of 12 mm was used in the experiments. This diameter was selected to ensure a balance between processing speed and precision. The use of a ball end mill in this research is based on its wide application in the 3D machining of complex geometrical surfaces, such as orthoses and molds, due to the tool’s ability to precisely process curved surfaces that correspond to those specified by the CAD model. Additionally, it facilitates efficient material removal without causing vibrations or compromising surface finish quality. By modifying the proposed algorithm, it is also possible to implement a flat-bottom tool, particularly for rough machining processes. However, the drawback of using a flat-bottom tool in 3D surface machining is the need for its rotation, which is very complex and prone to creating burrs. Therefore, this type of tool, as well as tools with more complex geometries, were not analyzed in this study, as they are rarely used in the production of orthoses and molds for orthoses.

The lead angle of the ball end mill can affect the generation of toolpaths [75], as changing the lead angle can improve processing efficiency in certain situations, especially when specific geometries or materials are present that require a customized machining approach. However, the proposed approach uses normal vectors as the basis for defining the tool direction. This means that the tool is always oriented perpendicular to the surface, eliminating the need for additional lead angle adjustments. As a result, the lead angle is effectively 0. This approach simplifies the management of the toolpath without the need for additional calculations arising from changes in the lead angle.

The processing parameters were as follows:

• Cutting speed: 3000 mm/min;
• Cutting depth (finishing): 2 mm;
• Spindle speed: 16,000 RPM.

The specified parameters were used in climb milling to achieve the best possible surface quality. In contrast to conventional milling, climb milling yields better results when machining softer materials, as it generates less heat on the surface during processing. This was particularly important in this study, which involved the use of soft polymer materials.

3.2.3. Results

Molds for orthopedic insoles were made from expanded polystyrene (XPS), while molds for spinal and ankle–foot orthoses were produced from expanded polyurethane (PU) (Figure 12 and Figure 13).

Although polyurethane and polystyrene are relatively soft materials that allow for faster processing, careful control of the parameters is essential to prevent overheating or surface damage. Therefore, the machining parameters varied empirically during initial testing depending on the characteristics of the two types of materials used, as well as the shape and volume of the workpieces. In the end, the previously stated final parameter values were taken as representative because they achieved the desired surface quality and processing efficiency for both types of materials. However, these parameters can still be considered initial and will be additionally refined and optimized in future studies.

The selected molds were those for which the simulation process confirmed that machining could be performed with the tool oriented according to the normal vector of the surface, given the technical specifications of the machine tool. In situations where there was a need to reduce the volume of the workpiece before the finishing process, rough machining was applied (left side of Figure 12). In this research, two methods were used for rough machining. In both methods, the tool enters the material vertically downward, parallel to the Z-axis. When generating the tool paths, the process begins by positioning the tool on the outermost (furthest from the axis of rotation) surface of the workpiece. The tool then moves inward until it touches the first layer of the positioned model, which, in this case, is represented by voxels. Each layer of material is defined by a certain height.

The first method is based on the principle of removing material layer by layer, with the workpiece rotating around the C- (Y-) axis. This process begins with the formation of a virtual cylinder around the axis of rotation, whose radius corresponds to the maximum distance between the workpiece and the voxelized model. A characteristic of this method is determining the distance from a point on the edge of the represented cylinder to the surface of the model that the tool passes through at a given position. This distance is divided by the layer height to calculate the number of passes required to remove the material. The depth of each tool pass is adjusted to the material of the workpiece, as well as to the processing conditions (type of tool and machining parameters). The material removal process begins along the axis of rotation, with the tool moving parallel to the axis of rotation from the minimum Y-axis to the maximum Y-axis coordinate of the workpiece. After completing the pass along this axis, the workpiece is rotated by a predefined angle around the C-axis, and the process is repeated. After each pass, the tool is lowered by the layer height towards the axis of rotation (C-axis), thus starting a new cycle of processing the next layer of material. This process is repeated until all layers are processed. In the end, the result is a roughly shaped model that follows the general dimensions of the CAD model but still contains a minimal amount of excess material that will be removed in the finishing process. This approach is efficient because it allows continuous rotation of the workpiece, thereby enabling material removal from all directions, resulting in more uniform removal of excess material.

The second method is based on generating a set of rotation angles for the A- and C-axes. These rotation angles depend on the complexity of the model’s geometry. In this research, nine combinations of workpiece rotation were used—three rotations around the C-axis (each by 120°) and three angles around the A-axis (−30°, 0°, +30°). Each combination of rotation angles is used for material removal by the ZIG-ZAG method, where for each set rotation angle, a top-view projection of the workpiece on the XY plane of the machine is generated. This projection enables precise determination of the region to be machined, defining the minimum and maximum coordinates on the X- and Y-axis. For this region, the layers of material to be removed are then defined. A characteristic of this method is that for a given starting point, a plane with a normal vector parallel to the Z-axis is determined, and the distance is calculated as the orthogonal distance from the underlying voxel to that plane. As in the first method, this distance is then divided by the layer height to calculate the number of passes needed for material removal. The material removal process is repeated for each combination of rotation angles until all excess material is removed. Both rough machining methods enable efficient and uniform removal of large amounts of material from all directions.

The primary criterion for evaluating the quality of the proposed toolpath generation method was the accuracy of physically replicating the original CAD model. The approach for digitizing the produced models and comparing them with the original CAD models was chosen. It included extracting the produced mold model from the surrounding material and placing it on a rotary table with a calibration object for scanning. Precise comparison of the physical models with the CAD models was ensured using the ATOS X5^® (GOM Gmbh—ZEISS Group, Braunschweig, Germany) 3D scanning system (Figure 14). Based on the surfaces obtained from multiple views, a single model was formed using the positioned calibration object.

The CAD model and the triangle mesh generated through digitization were imported into CloudCompare^® (ver. 2.6.1) software. The Iterative Closest Point (ICP) method was applied to register and minimize deviations between the corresponding surfaces. In the observed test models of orthopedic insoles, an average deviation of 0.47 mm across all compared CAD models and scanned surfaces was determined. One such example is presented in Figure 15, where a mean deviation of 0.26 mm between the physical and CAD models was observed.

The deviation of 0.47 mm is acceptable and falls within the commonly accepted criteria in practice, where deviations for orthopedic insoles should be less than 2 mm [76]. This was also the reason why the insole mold model was used for verifying the results in this work, as stricter criteria are applied to them in practice compared to, for example, spinal orthoses, where deviations of up to 10 mm are allowed. The results obtained were in line with expectations, confirming the validity of the method and providing a solid foundation for further improvements.

4. Conclusions

In this work, a new innovative method for the automatic generation of toolpaths for machining complex geometric surfaces on a five-axis machine has been proposed, which can be applied in the orthotics and prosthetics industry. The method is implemented by combining various algorithms, beginning with the discretization of the model into a voxel grid that defines its structure. Based on the generated grid, sequences of points are formed, from which equally spaced contours are then generated. Through a specially developed approach, the realized contours are filtered and organized to achieve the desired toolpaths.

The proposed method has proven to be stable, robust, and reliable in calculating toolpaths, ensuring consistent quality in the fabrication of orthosis molds. This stability is crucial, as irregularities in the toolpath can lead to inadequate orthosis production, negatively affecting the safety, comfort, and functionality of the device.

One of its main advantages is the ability to quickly and efficiently generate toolpaths for complex shapes. Traditional CAM software often requires detailed customization for each individual surface of the model, whereas the proposed algorithm automatically generates toolpaths based on the voxelization of the 3D model and wave propagation through the voxel grid, significantly reducing the time needed for toolpath preparation and manual processing of the workpiece models. The proposed approach allows for flexibility in choosing the voxel size, enabling a balance between precision and processing time, which are inversely proportional. Smaller voxels increase the accuracy of the toolpaths, while larger voxels reduce processing time. This adaptability is particularly beneficial in the orthotics industry, where the machining process must be quickly adjusted to accommodate different shapes and sizes of orthoses.

Although the method has been successfully tested on orthosis mold models with simpler geometries, further research is needed to confirm its applicability to more complex shapes and materials. Specifically, we aim to explore its application to other types of products that require similar surface precision and complexity. Future research will focus on several areas, as presented hereinafter.

First, we will analyze processing time and surface roughness (Ra) parameters in greater detail, as they were only roughly considered at this stage. This analysis aims to optimize the algorithm, particularly the toolpaths, for different materials and geometries, enhancing the quality of results while further reducing computation and processing time without compromising quality.

Second, to assess the algorithm’s adaptability across a wider range of industrial applications, we will extend the testing to more complex 3D orthosis models, various types of machines and materials, and explore the possibility of integrating the proposed method with other technologies, such as 3D scanning and additive manufacturing.

Third, we will conduct further comparisons with commercial CAM solutions to provide quantitative data on the advantages of the proposed algorithm in terms of precision, processing time, and surface quality.