An Efficient and Adaptable Path Planning Algorithm for Automated Fiber Placement Based on Meshing and Multi Guidelines

Path planning algorithms for automated fiber placement are used to determine the directions of the fiber paths and the start and end positions on the mold surfaces. The quality of the fiber paths determines largely the efficiency and quality of the automated fiber placement process. The presented work investigated an efficient path planning algorithm based on surface meshing. In addition, an update method of the datum direction vector via a guide-line update strategy was proposed to make the path planning algorithm applicable for complex surfaces. Finally, accuracy analysis was performed on the proposed algorithm and it can be adopted as the reference for the triangulation parameter selection for the path planning algorithm.


Introduction
Fiber-reinforced polymers (FRPs), especially carbon fiber reinforced polymers (CFRPs), are widely used in the aerospace and automobile industry as well as other fields because of their high specific strength, high specific modulus, excellent corrosion-and fatigue-resistance, and outstanding designability [1]. Automated fiber placement combines special manufacturing equipment with computerized numerical control (CNC) systems, which enables good forming quality, strong adaptability, and high efficiency. As an advanced composite manufacturing technology of large and complex components, automated fiber placement has become one of the most advanced and cutting-edge technologies for composite forming [2]. Path planning algorithms for automated fiber placement are used to specify the direction of a fiber path and determine the start and end positions of the fiber on the mold surface. The quality of the fiber path plays a crucial role for the efficiency and quality of the fiber placement process.
At present, the typical fiber path planning algorithms include the geodesic method and the meshing method based on parametric surfaces. The path planning algorithm of geodesic method based on parametric surfaces includes analytical and numerical methods [3]. The fiber path, which is obtained via the numerical solution of geodesics, is more aligned with practical engineering requirements. Lewis et al. [4] proposed a path planning method based on natural path, which is also widely used for path planning of tape laying. Shirinzadeh et al. [5][6][7] put forward the method of intersecting lines and surfaces to build the initial path, which improves the adaptability of the path to the surface curvature. However, there are certain disadvantages with the numerical solution of geodesics, e.g., the difficulty of finding the solution and the low efficiency associated with the calculation [8,9]. The meshing method, which uses the polygonal meshes to describe the complex parametric surfaces,

Efficiency Analysis of the Proposed Algorithm
A traditional path planning algorithm based on the geodesic method of parametric surface generates fiber paths by solving the direction of each point on the surface. It mainly uses the two geometric numerical solution operations, i.e., the intersection between a surface and a plane and the parallel offset of a curve. Essentially, the parallel offset of a curve is realized by the intersection of the surface at the sampling points and the corresponding offset direction planes.
Commercial CAD / CAM softwares generally use NURBS (Non-Uniform Rational B-Splines) for modelling, and the NURBS surface is represented by the following parameter equations: where Ω is the domain of s, , is the control point, , is the weight, , ( ) is the basis function of the pth-order B-spline in u direction, and , ( ) is the basis function of the qth-order B-spline in v direction.
Let the plane be represented by an implicit surface equation: where Ω is the domain of h.
Combining Equation (1) Intersection line l is a plane curve in the plane of parameter fields for u and v. If the surface s is the p × q-th order, the intersection line l is the p × q-th order. Because of the low efficiency of high-

Efficiency Analysis of the Proposed Algorithm
A traditional path planning algorithm based on the geodesic method of parametric surface generates fiber paths by solving the direction of each point on the surface. It mainly uses the two geometric numerical solution operations, i.e., the intersection between a surface and a plane and the parallel offset of a curve. Essentially, the parallel offset of a curve is realized by the intersection of the surface at the sampling points and the corresponding offset direction planes.
Commercial CAD/CAM softwares generally use NURBS (Non-Uniform Rational B-Splines) for modelling, and the NURBS surface is represented by the following parameter equations: where Ω s is the domain of s, P i,j is the control point, W i,j is the weight, N i,p (u) is the basis function of the pth-order B-spline in u direction, and N j,q (v) is the basis function of the qth-order B-spline in v direction.
Let the plane be represented by an implicit surface equation: where Ω h is the domain of h.
Combining Equation (1) with Equation (3) yields the intersection of plane h and surface s in domain Ω = Ω h ∩ Ω s . The parameters u and v of the intersection line satisfy the equation for the intersection line [14]: Intersection line l is a plane curve in the plane of parameter fields for u and v. If the surface s is the p × q-th order, the intersection line l is the p × q-th order. Because of the low efficiency of high-precision floating-point calculations, the frequent and high-precision solution of the intersection equation l consumes a lot of resources. This decreases both efficiency and stability of the numerical algorithms for the intersection between surfaces and planes. In the path planning algorithm based on the parametric surfaces, the intersection between surfaces and planes needs to be solved frequently. This decreases the efficiency of the path calculation, and for complex surfaces, the equation order is too high to be solved, which leads to the failure of path planning.
For triangular mesh surfaces, the surface is approximately a plane within the mesh cell [15]. The intersection between the surface and the plane can be converted into an intersection between planes, and the straight line in the plane is the result of the intersection.
The intersection line equation of a triangular patch A i (x, y, z) = 0 and the plane h is: The intersection line l is the first order equation about t, and the solution complexity decreases significantly. Since the triangulation algorithm will generate discrete grid planes of different sizes with different model complexity, the algorithm will not increase the order and complexity of equation solution during path generation due to the increase of model complexity. Therefore, the fiber path planning algorithm based on triangular mesh surface can obtain stable numerical solution quickly and efficiently.

Triangulation Algorithm for Parametric Surface
In the triangulation algorithm used in this paper, the edge and inner surfaces of the cell surface are approximated by straight line segments. Hence, the surface is discretized into area strips, and it is further divided into plane triangles. The triangulation results for the surface are generally given as strips of triangles. A strip of triangles is a list of points, such that any three consecutive points define a triangle. A parametric surface can be approximated by a series of triangle strips.
Due to the approximation process of replacing a curved surface with plane surfaces, the discrete error occurs in the triangulation process for parametric surface. This error is mainly reflected in the distance between the meshed plane surfaces and the original parametric surface. Two triangulation parameters are used to constrain the approximation error: (a) D l2s : The maximum distance from the straight-line segment of the discrete triangular area to the original parametric surface. (b) L seg : The maximum length of the straight-line segment in the discrete triangular area.

Sub-surface Boundary Splicing and Surface Topology Reconstruction
NURBS curves and surfaces, which are widely used in CAD/CAM software, have exact mathematical expressions, strong expression ability, good quality, and are easy to control. However, for complex surfaces, NURBS surfaces need to be split, spliced and trimmed. A complex surface is usually composed of multiple cellular surfaces. To reconstruct the topology of the whole complex surface, in addition to the mesh reconstruction of the cellular parameter surfaces, boundary splicing between cellular patches should be performed to obtain the global geometric topology.

Topology Reconstruction Algorithm for the Triangular Mesh of a Cellular Patch
As mentioned above, a cellular surface can be split into feature sampling points. To restore the surface information, it is necessary to obtain the segmentation results and establish the data structure for the triangular meshes by the vertex aggregation algorithm for triangulation (Algorithm 1). According to the right-hand rule, the vertices in the current discrete cell are stored to form the Point ID List of the current cell parameter surface. 3: Loop all points in the Point ID List 4: Every three points in the point table form a triangle patch to summarize the Face ID list and store the indexes of the three points of the current triangle patch.

5:
Point ID List update, add triangle patch ID index. 6: Build the edge ID list, update the two-point indexes of the edge, update the Point ID list to add the edge index, update the Face ID list to add the included edge index. 7: Calculate the normal vector of the current triangular patch according to the right-hand rule, and update it in the Face ID list.
Through the above process, the local Point ID list, Edge ID list and Face ID list are established, which reconstruct the topology information for the current NURBS cellular surface.

Algorithm for Subsurface-Boundary Splicing
Because the mesh discretization results for each cellular parameter surface are independent of each other, and two adjacent cellular surfaces share the same edge, there are duplicate vertices for the adjacent cellular surfaces. It is necessary to remove the duplicate vertices. The subsurface-boundary splicing algorithm includes removing duplicate vertices and updating the global index of the vertices in all the cellular surfaces.
The brute force algorithm, which traverses the whole Point ID List for the duplicate vertices with the same coordinates with a given vertex and hence delete the duplicate ones, is the most straightforward method. Suppose that the surface is made up of k quadrilateral cellular NURBS surfaces with similar area and all the NURBS surfaces are smooth. Following the triangulation, it was assumed that the triangulation results are all given as strips of triangles to obtain a uniform triangular network. If the number of vertices on the boundary is a and b, the number of vertices on the cellular surface is a × b, and the number of all vertices is N = k × a × b. The time complexity of the brute force algorithm is O(N 2 ).
However, all the duplicate vertices are located on the boundary line because they are formed by two adjacent cellular surfaces sharing a common edge. For the vertices on the boundary line, they are in a semi-closed state and not surrounded by all triangles. On the other hand, the vertices inside the surfaces are in a fully-closed state, surrounded by several triangles (Figure 2). When the vertex is in the fully-closed state, the number of adjacent patches is equal to the number of edges, while in the semi-closed state, the number of patches is not equal to the number of edges. Therefore, all triangle vertices can be divided into two types: boundary semi-closed vertices and internal fully-closed vertices. Thanks to this feature, the boundary vertex set can be filtered out. Then, duplicate vertices can be removed from the boundary vertex set, and Point ID list, Face ID list and Edge ID list can be updated simultaneously. The time complexity of the algorithm is O (k × (a + b) 2 ), which is far less than O (N 2 ) of the brute force algorithm, and the efficiency of the algorithm is greatly improved. At the end of the above process, both sub-surface boundary splicing and surface topological relation reconstruction are completed. This enables fast search of adjacent triangles and efficient path calculation. The topological information of the discrete mesh surfaces is contained in the Point ID list, Face ID list and Edge ID list. The forms and requirements of the three data structure list are shown in Table 1 [13]. Given the number of a triangle patch, one can quickly find the global index and the global index of the three edges of the three vertices on the triangle patch. De-duplication algorithm based on the vertex closure checking was described in Algorithm 2. The time complexity of the algorithm is O (k × (a + b) 2 ), which is far less than O (N 2 ) of the brute force algorithm, and the efficiency of the algorithm is greatly improved.
At the end of the above process, both sub-surface boundary splicing and surface topological relation reconstruction are completed. This enables fast search of adjacent triangles and efficient path calculation. The topological information of the discrete mesh surfaces is contained in the Point ID list, Face ID list and Edge ID list. The forms and requirements of the three data structure list are shown in Table 1 [13]. Given the number of a triangle patch, one can quickly find the global index and the global index of the three edges of the three vertices on the triangle patch.

Main Path Planning Algorithm Based on Guidelines on Triangular Mesh Surfaces
The path planning process needs to generate the main motion paths of the head of the automated fiber placement equipment and the corresponding fiber paths which are offset by the main motion Materials 2020, 13, 4209 6 of 18 paths. The fiber path generation is relatively simple and this paper focuses on the main path planning algorithm. The guidelines are extracted from the CAD model of the mold of the to-be-fabricated FRP component to reflect the skeleton and unique appearance of the component. According to the guide-lines, the direction vector of the main path can be calculated and adjusted adaptively, so that the main path and the corresponding fiber paths can comply with the shape of the component. It is beneficial to improve the mechanical performances of the FRP component while satisfying the constraints of the minimum steering radius of the fiber materials, even distribution of the cutting points, etc.

Main Motion Path Generation Algorithm
According to the fiber placement process, three geometric parameters (initial starting point P 0 , guide line L, parametric surface Π) as well as the ply angle θ should be provided as input to the main motion path generation algorithm, which outputs the main motion path l i . Basically, the algorithm needs to calculate the direction vectors and subsequently generate the continuous points on the main motion path, as described in Algorithms 3 and 4, respectively.
Algorithm 3: The direction vector calculation algorithm for the main motion path.
Step 1: Find the projection point P 0 ' for the initial starting point P 0 on the triangle patch A of the triangular mesh surfaces (triangulation of the parametric surface Π) and obtain the normal vector n of triangle patch A. Then, redefine P 0 ' as the starting point for the current main motion path.
Step 2: Calculate the tangent vector t' for the projection point P 0 " of P 0 ' on the guide line L, and generate the parallel vector t of t' at point P 0 '.
Step 3: Calculate the orthogonal vector k of the normal vector n and the parallel tangent vector t, where k = n × t.
Step 4: Calculate the orthogonal vector m of normal vector n and vector k, where m = k × n. The vector m is the datum direction vector for the current path point, based on the guide-line and corrected by the curvature feature of the surface, where the path point is located (as shown in Figure 3).
The path planning process needs to generate the main motion paths of the head of the automated fiber placement equipment and the corresponding fiber paths which are offset by the main motion paths. The fiber path generation is relatively simple and this paper focuses on the main path planning algorithm. The guidelines are extracted from the CAD model of the mold of the to-be-fabricated FRP component to reflect the skeleton and unique appearance of the component. According to the guidelines, the direction vector of the main path can be calculated and adjusted adaptively, so that the main path and the corresponding fiber paths can comply with the shape of the component. It is beneficial to improve the mechanical performances of the FRP component while satisfying the constraints of the minimum steering radius of the fiber materials, even distribution of the cutting points, etc.

Main Motion Path Generation Algorithm
According to the fiber placement process, three geometric parameters (initial starting point P0, guide line L, parametric surface Π) as well as the ply angle θ should be provided as input to the main motion path generation algorithm, which outputs the main motion path li. Basically, the algorithm needs to calculate the direction vectors and subsequently generate the continuous points on the main motion path, as described in Algorithms 3 and 4, respectively.
Algorithm 3: The direction vector calculation algorithm for the main motion path.
Step 1: Find the projection point P0' for the initial starting point P0 on the triangle patch A of the triangular mesh surfaces ∑ (triangulation of the parametric surface Π) and obtain the normal vector n of triangle patch A. Then, redefine P0' as the starting point for the current main motion path.
Step 2: Calculate the tangent vector t' for the projection point P0'' of P0' on the guide line L, and generate the parallel vector t of t' at point P0'.
Step 3: Calculate the orthogonal vector k of the normal vector n and the parallel tangent vector t, where k = n × t.
Step 4: Calculate the orthogonal vector m of normal vector n and vector k, where m = k × n. The vector m is the datum direction vector for the current path point, based on the guide-line and corrected by the curvature feature of the surface, where the path point is located (as shown in Figure  3).
Step 5: Using normal vector n as the rotation axis, rotate the datum direction vector m for θ degree to obtain the laying direction vector d for the current path point P0'.  Step 5: Using normal vector n as the rotation axis, rotate the datum direction vector m for θ degree to obtain the laying direction vector d for the current path point P 0 '.
In the current triangle patch A, use point P 0 and vector d to construct a straight line and hence obtain the intersection point P 1 of the straight line and the three sides of triangle patch A. The intersection point P 1 is the next path point, and all the path points can be obtained after continuous iteration for the main motion path generation.
Algorithm 4: The continuous path point generation algorithm.
Step 1: Take the starting point P i and direction vector d i in the current triangle patch A i , construct a ray P a (λ), which is defined as the solution line for the path point. Its parameter equation is: Step 2: The vertices V a and V b of the triangle patch form a straight line P b (λ), and the parameter equation is: Step 3: Using P a (λ) = P b (λ), we can find the unique solution λ: If 0 ≤ λ ≤ 1 is true, this means that the next path point is on the current sideline. If it is not on the current sideline, then take two points to form a straight line for the calculation. Next, take V c V b and V a V c to form a straight line for the calculation. Repeat steps 1, 2, and 3 to get the next path point.
The solution line for the path point may overlap with the three side lines of the triangle in step 3, and the next path point P i+1 cannot be calculated using the above steps. At this time, the endpoint of the edge V m V k (m, k ∈ {a, b, c}) of the triangle, where the point P i is located, is used as the next path point P i+1 . The endpoint selection rules are described in Equation (11), where d is the laying direction vector.
Step 4: According to the topological relationship for the triangular mesh surface, obtain the adjacent triangles on the other side of the edge, where P i+1 is located, and update the patch index as A i+1 .
Step 5: Update the normal vector n i+1 . When the next path point P i+1 lies on the triangle edge or vertex, use the area weighing method to reduce the deviation.
In Equation (12), m is the total number of adjacent triangles to which P i+1 belongs, and A k is the area of the triangles. The area can be determined using A k = AB×AC 2 , where points A, B, C are the three vertices of a triangle A k .
Step 6: Update the parallel tangent vector t i+1 , and solve the direction vector d i+1 of the next path point P i+1 using Algorithm 3. The projection vector d i+1 of d i+1 on patch A k+1 is the laying-direction vector of point P i+1 .
Step 7: Repeat Step 1 to Step 6 to update and calculate the next path point until reach the triangular mesh surface boundary and no adjacent patch can be found and the partial path on the one single direction is constructed.
Step 8: Go back to the path initial point P 0 , take the inverse of the laying-direction vector d 0 and repeat Step 1 to Step 7 until the whole main motion path is completely constructed. The continuous path point generation algorithm is shown in Figure 4. After completing the above steps, the main motion path can be derived using the cubic spline interpolation algorithm.

Guide-line Update Algorithm for Complex Surfaces
If the curvature of the surface of certain components changes substantially, the path, which was calculated and planned using a single guide line, cannot meet the requirements for the laying ability in each surface area. This leads to wrinkling of the fiber during the laying process, which degrades the mechanical properties of the final fabricated components. To address this problem, this paper proposed a guide-line algorithm for the update of the tangent vector t and datum direction vector m for the path points according to the surface shape and the distribution of multi guide-lines, as demonstrated in Figure 5 and described in Algorithm 5.

Guide-line Update Algorithm for Complex Surfaces
If the curvature of the surface of certain components changes substantially, the path, which was calculated and planned using a single guide line, cannot meet the requirements for the laying ability in each surface area. This leads to wrinkling of the fiber during the laying process, which degrades the mechanical properties of the final fabricated components. To address this problem, this paper proposed a guide-line algorithm for the update of the tangent vector t and datum direction vector m for the path points according to the surface shape and the distribution of multi guide-lines, as demonstrated in Figure 5 and described in Algorithm 5.

Guide-line Update Algorithm for Complex Surfaces
If the curvature of the surface of certain components changes substantially, the path, which was calculated and planned using a single guide line, cannot meet the requirements for the laying ability in each surface area. This leads to wrinkling of the fiber during the laying process, which degrades the mechanical properties of the final fabricated components. To address this problem, this paper proposed a guide-line algorithm for the update of the tangent vector t and datum direction vector m for the path points according to the surface shape and the distribution of multi guide-lines, as demonstrated in Figure 5 and described in Algorithm 5.  Step 1: Project the current path point P i on each guide-line L i to obtain projection points A, B, C and D. Connect the projection points to form a polygonal section plane Ω, the four sides of the polygon on the section plane Ω are a, b, c and d.
Step 2: Project the current path point P i to the edge of the polygon, and find the projection points P 1 i , P 2 i , P 3 i and P 4 i .
Step 3: If the projection point is on the extension line of the edge, the edge is excluded. As shown in Figure 5b, the projection point P 4 i is not on the edge d, so the edge d is not considered in the algorithm below.
Step 4: Calculate the distance from P i to the projection points on the not excluded edge lines, where dist = P i P x i .
Step 5: Determine the minimum distance dist and its edge x. The two guide-lines L a and L b , at the end of x, are the guide-lines for the current path point P i .
Step 6: As shown in Figure 5c and Equation (14), calculate the tangent vectors on L a and L b , respectively. The smooth transition for the path direction vector from L a to L b is realized using the distance-weighing method: where dist 1 and dist 2 are the distances of P i to L a and L b . Update the tangent vector t and calculate the datum direction vector m of P i in Algorithms 3 and 4. The path, which is generated by Algorithm 5, is more adaptive to the changing curvature of the surface and the scalability is improved. Figure 6 shows the paths generated based on single guideline and multi-guidelines for a panel and a curved surface models. It can be seen that the fiber direction transits smoothly between the guide lines and this makes the placed fibers adapt to the shape of the moulds of the final parts.
Materials 2020, 13, x FOR PEER REVIEW 9 of 19 Algorithm 5: Guideline update algorithm for complex surfaces.
Step 1: Project the current path point Pi on each guide-line Li to obtain projection points A, B, C and D. Connect the projection points to form a polygonal section plane Ω, the four sides of the polygon on the section plane Ω are a, b, c and d.
Step 2: Project the current path point Pi to the edge of the polygon, and find the projection points , , and .
Step 3: If the projection point is on the extension line of the edge, the edge is excluded. As shown in Figure 5b, the projection point is not on the edge d, so the edge d is not considered in the algorithm below.
Step 4: Calculate the distance from Pi to the projection points on the not excluded edge lines, where = ‖ ‖.
Step 5: Determine the minimum distance dist and its edge x. The two guide-lines La and Lb, at the end of x, are the guide-lines for the current path point Pi.
Step 6: As shown in Figure 5c and Equation (14), calculate the tangent vectors on La and Lb, respectively. The smooth transition for the path direction vector from La to Lb is realized using the distance-weighing method: where dist1 and dist2 are the distances of Pi to La and Lb. Update the tangent vector t and calculate the datum direction vector m of Pi in Algorithms 3 and 4. The path, which is generated by Algorithm 5, is more adaptive to the changing curvature of the surface and the scalability is improved. Figure 6 shows the paths generated based on single guideline and multi-guidelines for a panel and a curved surface models. It can be seen that the fiber direction transits smoothly between the guide lines and this makes the placed fibers adapt to the shape of the moulds of the final parts.

Accuracy Analysis of the Generated Path Based on Surface Meshing
In the triangulation process in Section 2, D l2s and L seg are the main parameters that affect both the mesh density and the approximation accuracy. When the two parameters are set to larger values, the mesh density is small, the surface approximation accuracy is low, the subsequent path planning process data volume is small, and the algorithm efficiency is high. If the two parameter values are continuously reduced, on the other hand, the efficiency is lower.
It is generally believed that the path generated using the geodesic method on the parametric surface is used as the standard path when the planning accuracy and error of the algorithm need to be verified. The points on the discrete surface of the mesh are used to replace the path points on the original parametric surface. Furthermore, the normal vector of the path points on the triangular mesh surface can be used to replace the normal vector on the original parametric surface. As a result, a cumulative error can occur during the iteration process of the proposed path planning algorithm, which can cause the angle deviation from the design datum and the distance deviation from the path generated on the original parametric surface.
A complex surface with positive and negative curvature was adopted to conduct accuracy analysis of the proposed path planning algorithm, as shown in Figure 7. The surface consists of 29 independent cellular patches. Each cellular parameter surface represents a NURBS surface, and the order of each cellular parameter surface in both U and V directions is 6 degrees.

Accuracy Analysis of the Generated Path Based on Surface Meshing
In the triangulation process in Section 2, Dl2s and Lseg are the main parameters that affect both the mesh density and the approximation accuracy. When the two parameters are set to larger values, the mesh density is small, the surface approximation accuracy is low, the subsequent path planning process data volume is small, and the algorithm efficiency is high. If the two parameter values are continuously reduced, on the other hand, the efficiency is lower.
It is generally believed that the path generated using the geodesic method on the parametric surface is used as the standard path when the planning accuracy and error of the algorithm need to be verified. The points on the discrete surface of the mesh are used to replace the path points on the original parametric surface. Furthermore, the normal vector of the path points on the triangular mesh surface can be used to replace the normal vector on the original parametric surface. As a result, a cumulative error can occur during the iteration process of the proposed path planning algorithm, which can cause the angle deviation from the design datum and the distance deviation from the path generated on the original parametric surface.
A complex surface with positive and negative curvature was adopted to conduct accuracy analysis of the proposed path planning algorithm, as shown in Figure 7. The surface consists of 29 independent cellular patches. Each cellular parameter surface represents a NURBS surface, and the order of each cellular parameter surface in both U and V directions is 6 degrees.

Distance Deviation Analysis
The distance deviation of the path generated on the triangular mesh surface can be decomposed into the normal distance deviation and the geodesic distance deviation. In this paper, a uniform orthogonal test was carried out for the appropriate ranges of Dl2s and Lseg. The path was uniformly sampled (with a distance of 5 mm) to evaluate the distance deviation, as shown in Figure 8.
During the triangulation process, both Dl2s and Lseg did not reach 0, which means that the parameter value (without approximation error) could not be determined. Therefore, in the actual application, the Dl2s range was 0.2-1.0 mm, and the Lseg range was 20-200 mm. The experiment was carried out by the (5 ) orthogonal design. The path-generation time, the distance deviation, the normal distance deviation and the geodesic distance deviation were recorded.

Distance Deviation Analysis
The distance deviation of the path generated on the triangular mesh surface can be decomposed into the normal distance deviation and the geodesic distance deviation. In this paper, a uniform orthogonal test was carried out for the appropriate ranges of D l2s and L seg . The path was uniformly sampled (with a distance of 5 mm) to evaluate the distance deviation, as shown in Figure 8.
During the triangulation process, both D l2s and L seg did not reach 0, which means that the parameter value (without approximation error) could not be determined. Therefore, in the actual application, the D l2s range was 0.2-1.0 mm, and the L seg range was 20-200 mm. The experiment was carried out by the L 25 5 6 orthogonal design. The path-generation time, the distance deviation, the normal distance deviation and the geodesic distance deviation were recorded.
The Euclidean distance deviation d from the sample point on the generated path to the standard reference path was calculated and decomposed into the normal distance deviation d N and the geodesic distance deviation d T . The Euclidean distance deviation d from the sample point on the generated path to the standard reference path was calculated and decomposed into the normal distance deviation dN and the geodesic distance deviation dT.
The project vector from the sample point to standard reference path is V, and the normal vector for sample point Pi on the original parametric surface is n. dN and dT can be calculated as follows: According to the The initial path point was located near the 350th sample point. According to Figure 9, the closer the sample point was to the initial point, the smaller was the distance deviation. During path generation, the normal vector for the path points on the triangular mesh surface was used to (approximately) replace the normal vector of the path points on the original parameter surface. Therefore, the direction datum of the path point was deviated, and it caused geodesic distance deviation between the generated path and the standard reference path. It also shows that the geodesic distance deviation dT was the main deviation and it increased as the parameters Dl2s and Lseg increased.
To analyze the relationship between the parameters (i.e., Dl2s and Lseg) and the distance deviation of the generated path and the algorithm efficiency, the mean distance deviation, the normal mean distance deviation, the geodesic mean distance deviation, the maximum distance deviation, the maximum normal distance deviation, the maximum geodesic distance deviation, and the path generation time were calculated. This was done when the generation time started from the beginning of the surface triangulation to the end of the path generation, as shown in Table 2. The project vector from the sample point to standard reference path is V, and the normal vector for sample point P i on the original parametric surface is n. d N and d T can be calculated as follows: According to the L 25 5 6 orthogonal design, 25 sets of experiments were carried out. Four of them, which were D l2s = 0.2 and L seg = 65, D l2s = 0.4 and L seg = 110, D l2s = 0.6 and L seg = 155 and D l2s = 0.8 and L seg = 200, were selected to analyze the distribution of distance deviation along the generated paths.
The initial path point was located near the 350th sample point. According to Figure 9, the closer the sample point was to the initial point, the smaller was the distance deviation. During path generation, the normal vector for the path points on the triangular mesh surface was used to (approximately) replace the normal vector of the path points on the original parameter surface. Therefore, the direction datum of the path point was deviated, and it caused geodesic distance deviation between the generated path and the standard reference path. It also shows that the geodesic distance deviation d T was the main deviation and it increased as the parameters D l2s and L seg increased.
To analyze the relationship between the parameters (i.e., D l2s and L seg ) and the distance deviation of the generated path and the algorithm efficiency, the mean distance deviation, the normal mean distance deviation, the geodesic mean distance deviation, the maximum distance deviation, the maximum normal distance deviation, the maximum geodesic distance deviation, and the path generation time were calculated. This was done when the generation time started from the beginning of the surface triangulation to the end of the path generation, as shown in Table 2.    The distribution of mean distance deviation are shown in Figure 10, while the maximum distance deviation are shown in Figure 11, and the generation time are shown in Figure 12.
Supplementary Video S1 further demonstrates the high efficiency of the proposed algorithm, which can complete the path planning for one layer of a complex surface in just only tens of seconds. The distribution of mean distance deviation are shown in Figure 10, while the maximum distance deviation are shown in Figure 11, and the generation time are shown in Figure 12.    The distribution of mean distance deviation are shown in Figure 10, while the maximum distance deviation are shown in Figure 11, and the generation time are shown in Figure 12.   According to Figures 10 and 11, for different parameters, the distance deviation mainly consisted of geodesic distance deviation, and the effect of normal distance deviation on the overall distance deviation is relatively small. Meanwhile, it can be observed that D l2s and L seg restricted each other in the triangulation process. When the two parameters cannot be satisfied simultaneously, the algorithm will adopt the parameter that makes the mesh more precise. For example, when L seg is 20 mm, the triangulation algorithm generated the same triangular mesh surfaces for a D l2s of 0.4 mm, 0.6 mm, 0.8 mm and 1.0 mm as for the D l2s of 0.2 mm. Table 2 shows that, when D l2s exceeds 0.8mm and L seg exceeds 155 mm, the mean distance deviation surpasses 1mm, and the maximum distance deviation exceeds 2 mm. When D l2s is between 0.2 and 0.6mm, and L seg is 65 to 155 mm, the average distance deviation remains within 1mm, while the maximum distance deviation is within 2 mm. The generation time of the algorithm is within 0.5 s. According to Figures 10 and 11, for different parameters, the distance deviation mainly consisted of geodesic distance deviation, and the effect of normal distance deviation on the overall distance deviation is relatively small. Meanwhile, it can be observed that Dl2s and Lseg restricted each other in the triangulation process. When the two parameters cannot be satisfied simultaneously, the algorithm will adopt the parameter that makes the mesh more precise. For example, when Lseg is 20 mm, the triangulation algorithm generated the same triangular mesh surfaces for a Dl2s of 0.4 mm, 0.6 mm, 0.8 mm and 1.0 mm as for the Dl2s of 0.2 mm. Table 2 shows that, when Dl2s exceeds 0.8mm and Lseg exceeds 155 mm, the mean distance deviation surpasses 1mm, and the maximum distance deviation exceeds 2 mm. When Dl2s is between 0.2 and 0.6mm, and Lseg is 65 to 155 mm, the average distance deviation remains within 1mm, while the maximum distance deviation is within 2 mm. The generation time of the algorithm is within 0.5 s.

Angle Deviation Analysis
To ensure the orthogonal arrangement of fibers in different layers, the angle distribution between different layers is strictly defined. This ensures quasi-isotropic mechanical behavior for the components [3]. To analyze the angle deviation between the generated paths and the design datum, an angle deviation analysis was performed in this section.
The generated path was sampled uniformly with a step of 5 mm. For each sample point, the forward direction of the current path and the design direction datum of the path were calculated using both the guide-line and surface normal vector at the point. Subsequently, the angle deviation was obtained. The forward direction vector at each sample point is the tangent vector of the sample point on the path. The datum direction vector can be calculated using:

Angle Deviation Analysis
To ensure the orthogonal arrangement of fibers in different layers, the angle distribution between different layers is strictly defined. This ensures quasi-isotropic mechanical behavior for the components [3]. To analyze the angle deviation between the generated paths and the design datum, an angle deviation analysis was performed in this section.
The generated path was sampled uniformly with a step of 5 mm. For each sample point, the forward direction of the current path and the design direction datum of the path were calculated using both the guide-line and surface normal vector at the point. Subsequently, the angle deviation was obtained. The forward direction vector d i at each sample point is the tangent vector of the sample point on the path. The datum direction vector d i can be calculated using: where n i is the normal vector of the sample point on the surface, and t i is the tangent of the projection point for the sample point on the guide-line. The angle deviation δ i is calculated using: The four sets of experiments selected in Section 4.1 were also adopted here to analyze the distribution of the angle deviation.
As mentioned above, the initial path point was near the 350th sampling point. According to Figure 13, for different triangulation parameters, the angle deviation does not depend on the distance between the sampling point and the initial path point but depend on the curvature of the original surface. The larger the curvature of the surface, where the sample points are located on the path, the greater is the angle deviation. This confirms that the angle deviation is due to the approximation of the normal vector of the triangular mesh surface during the execution of the algorithm, and there is no cumulative error.
The four sets of experiments selected in Section 4.1 were also adopted here to analyze the distribution of the angle deviation.
As mentioned above, the initial path point was near the 350 th sampling point. According to Figure 13, for different triangulation parameters, the angle deviation does not depend on the distance between the sampling point and the initial path point but depend on the curvature of the original surface. The larger the curvature of the surface, where the sample points are located on the path, the greater is the angle deviation. This confirms that the angle deviation is due to the approximation of the normal vector of the triangular mesh surface during the execution of the algorithm, and there is no cumulative error. To study the effect of Dl2s and Lseg on the angle deviation of the path, the mean angle deviation, the mean square error, and the maximum angle deviation were calculated using the experimentally obtained data-see Table 3. To study the effect of D l2s and L seg on the angle deviation of the path, the mean angle deviation, the mean square error, and the maximum angle deviation were calculated using the experimentally obtained data-see Table 3.  The distribution of average angle deviation data is shown in Figure 14, and the maximum angle deviation data distribution is shown in Figure 15.  According to Figures 14 and 15, as the D l2s and L seg increased, both the mesh density and approximation accuracy of the triangular mesh surface decrease. In addition, the normal vector of the path points on the triangular mesh surface deviate significantly from the normal vector on the original parametric surface. Hence, the angle deviation increases. According to Table 3, when D l2s exceeds 0.6 mm and L seg exceeds 65 mm, the mean angle deviation surpasses 0.25 deg, and the maximum deviation is more than 1.4 deg.
Based on the above distance deviation analysis, when D l2s ranges between 0.2 and 0.6 mm and L seg is 65 to 110 mm, the mean distance deviation remains within 1mm. Furthermore, the maximum distance deviation stays within 2 mm, the mean angle deviation is less than 0.25 deg, and the maximum angle deviation is below 1.4 deg. At the same time, the path generation time stays within 0.5s. In this case, D l2s and L seg can be selected according to the complexity of the surface and the acceptable path error. Within the above range of the parameters, a high-precision triangular mesh surface and fiber path with small error can be obtained, while the generation efficiency of the algorithm is high.

Conclusions
To improve the efficiency of automated fiber path planning process, a new path planning algorithm based on meshing and multi guide-lines were investigated. The original parameter surface of the CAD model of the FRP component was discretized into triangular mesh surface via surface discretization and triangulation. Sub-surface boundary splicing and surface topology reconstruction algorithm was proposed, and both the computational complexity reduction and the efficiency improvement of the algorithm were analyzed. The proposed automated fiber path planning algorithm consists of a main motion path direction vector algorithm and a continuous path point generation algorithm. An updating method for the datum direction vector via the guide-lines update algorithm was also introduced for complex surfaces. It improves the laying ability of the fibers and surface adaptability for the planned path. Accuracy analysis was conducted to investigate the relationship between the triangulation parameters and distance deviation, angle deviation and algorithm efficiency. The analysis indicated that by choosing appropriate triangulation parameters, the fiber path can be generated with high accuracy and efficiency.
More research efforts in the future work should be devoted to conduct experiments to test the mechanical properties of the fabricated FRP components by using the multi-guide-line planned paths.