Path Planning for Laser Cladding Robot on Artiﬁcial Joint Surface Based on Topology Reconstruction

: Artiﬁcial joint surface coating is a hot issue in the interdisciplinary ﬁelds of manufacturing, materials and biomedicine. Due to the complex surface characteristics of artiﬁcial joints, there are some problems with e ﬃ ciency and precision in automatic cladding path planning for coating fabrication. In this study, a path planning method for a laser cladding robot for artiﬁcial joints surface was proposed. The key of this method was the topological reconstruction of the artiﬁcial joint surface. On the basis of the topological relation, a set of parallel planes were used to intersect the CAD model to generate a set of continuous, directed and equidistant surface transversals on the artiﬁcial joint surface. The arch height error method was used to extract robot interpolation points from surface transversal lines according to machining accuracy requirements. The coordinates and normal vectors of interpolation points were used to calculate the position and pose of the robot tool center point (TCP). To ensure that the laser beam was always perpendicular to the artiﬁcial joint surface, a novel laser cladding set-up with a robot was designed, of which the joint part clamped by a six-axis robot moved while the laser head was ﬁxed on the workbench. The proposed methodology was validated with the planned path on the surface of an artiﬁcial acetabular cup using simulation and experimentation via an industrial NACHI robot. The results indicated that the path planning method based on topological reconstruction was feasible and more e ﬃ cient than the traditional robot teaching method.


Introduction
Artificial joints are widely used in medical treatment, which can effectively help patients with joint injury to eliminate pain and rebuild joint functions [1,2]. With the aggravation of aging problems and the increase of joint injuries caused by accidents such as traffic accidents, the demand for artificial joints is increasing day by day [3]. Titanium alloys are commonly used to make artificial joint prosthesis, due to its superior mechanical properties, excellent corrosion resistance and good biocompatibility [4]. However, artificial joint prosthesis, as a body implant in direct contact with bone or tissue fluid, needs to have excellent biological activity [5], and the contact surface of the joint friction pair needs to have extremely high wear resistance [6]. In order to meet the special requirements of artificial joint application, titanium alloys are usually selected as the substrate, and one or more layers of high performance coating is fabricated on the surface of titanium alloys. In recent years, laser cladding as a novel and advanced surface modification technology has been gradually used to prepare ceramic coatings on the surface of titanium alloys to improve the wear resistance or impart bioactive properties [7,8].
However, the geometric shape of artificial joints is complex [9], so it is difficult to plan the laser cladding path. How to design a suitable laser cladding path on such a complex surface and improve the efficiency and quality of coatings preparation has become an urgent problem to be Figure 1a. The laser cladding process is shown in Figure 1b. The coaxial head remains stationary, while the experimental sample carried by the six-axis robot moves along the planned cladding path. The laser beam stays perpendicular to the surface of the experimental sample all the time.
Algorithms 2020, 13, x FOR PEER REVIEW 3 of 19 while the experimental sample carried by the six-axis robot moves along the planned cladding path. The laser beam stays perpendicular to the surface of the experimental sample all the time. The method proposed in this paper was implemented using processing software. Its architecture is shown in Figure 2. It mainly consisted of three modules: the CAD module, the path planning module and the simulation module. The CAD module offered two functions: graphic user interface for model handling and basic geometric operations for model transformation. The path planning module offered three functions: topological reconstruction of surface, cladding path planning and robot program generation. The simulation module was provided by the robot manufacturer.

The research methods
The overall procedure of laser cladding path planning proposed in this paper is shown in Figure  3. The goal was to generate an optimal laser cladding path plan for the freeform surface parts of artificial joints. The overall procedure included mainly three steps: Step 1: Analysis of artificial joint CAD model. The topological relation of the artificial joint surface was created by contacting feature vertices, edges and facets extracted from the artificial joint CAD model.
Step 2: Generation of the cladding paths according to the CAD model. The laser cladding path was automatically generated using the slicing algorithm based on a topological relation. The arch height error method was used to extract reasonable cladding points from a large number of intersections.
Step 3: Determination of the robot position and pose. The position and pose of the robot TCP were calculated according to the coordinates and normal vectors of the interpolation points. The The method proposed in this paper was implemented using processing software. Its architecture is shown in Figure 2. It mainly consisted of three modules: the CAD module, the path planning module and the simulation module. The CAD module offered two functions: graphic user interface for model handling and basic geometric operations for model transformation. The path planning module offered three functions: topological reconstruction of surface, cladding path planning and robot program generation. The simulation module was provided by the robot manufacturer.
Algorithms 2020, 13, x FOR PEER REVIEW 3 of 19 while the experimental sample carried by the six-axis robot moves along the planned cladding path. The laser beam stays perpendicular to the surface of the experimental sample all the time. The method proposed in this paper was implemented using processing software. Its architecture is shown in Figure 2. It mainly consisted of three modules: the CAD module, the path planning module and the simulation module. The CAD module offered two functions: graphic user interface for model handling and basic geometric operations for model transformation. The path planning module offered three functions: topological reconstruction of surface, cladding path planning and robot program generation. The simulation module was provided by the robot manufacturer.

The research methods
The overall procedure of laser cladding path planning proposed in this paper is shown in Figure  3. The goal was to generate an optimal laser cladding path plan for the freeform surface parts of artificial joints. The overall procedure included mainly three steps: Step 1: Analysis of artificial joint CAD model. The topological relation of the artificial joint surface was created by contacting feature vertices, edges and facets extracted from the artificial joint CAD model.
Step 2: Generation of the cladding paths according to the CAD model. The laser cladding path was automatically generated using the slicing algorithm based on a topological relation. The arch height error method was used to extract reasonable cladding points from a large number of intersections.
Step 3: Determination of the robot position and pose. The position and pose of the robot TCP were calculated according to the coordinates and normal vectors of the interpolation points. The

The Research Methods
The overall procedure of laser cladding path planning proposed in this paper is shown in Figure 3. The goal was to generate an optimal laser cladding path plan for the freeform surface parts of artificial joints. The overall procedure included mainly three steps: Step 1: Analysis of artificial joint CAD model. The topological relation of the artificial joint surface was created by contacting feature vertices, edges and facets extracted from the artificial joint CAD model.
Step 2: Generation of the cladding paths according to the CAD model. The laser cladding path was automatically generated using the slicing algorithm based on a topological relation. The arch height error method was used to extract reasonable cladding points from a large number of intersections.
Step 3: Determination of the robot position and pose. The position and pose of the robot TCP were calculated according to the coordinates and normal vectors of the interpolation points. The robot motion paths were obtained by connecting the interpolation points sequentially, which finally were converted into the robot program.
Algorithms 2020, 13, x FOR PEER REVIEW 4 of 19 robot motion paths were obtained by connecting the interpolation points sequentially, which finally were converted into the robot program.

Topological reconstruction of the artificial joint CAD model
A STL (Stereo Lithography) file, which is a CAD data exchange file with a simple format, is very suitable for simulating, designing and modifying artificial joint models [20]. 3D models of artificial hip joints were converted into STL files by the CAD module, as shown in Figure 4a. Figure 4b shows the triangular facets in the STL file, which were described using three sets of X, Y and Z coordinates for each vertex and a normal unit vector indicating the direction of the facet. The three vertices were arranged counterclockwise. It must be mentioned that the triangular facets in the STL file were isolated and had no relation to the surrounding triangles. In other words, the STL file contained only graphical location information, not topological information. Therefore, in order to extract geometric information more effectively in the CAD model for subsequent path planning, it was necessary to pre-process STL files and create the relations among points, edges and facets [21]. In order to create the overall topological information at one time, the topology access should meet the following basic rules: Rule 1: Given a triangular facet, get its three edges and three vertices. Rule 2: Given an edge, get its two adjacent triangular facets and two vertices. Rule 3: Given a vertex, get all triangular facets in which the vertex is located.

Topological Reconstruction of the Artificial Joint CAD Model
A STL (Stereo Lithography) file, which is a CAD data exchange file with a simple format, is very suitable for simulating, designing and modifying artificial joint models [20]. 3D models of artificial hip joints were converted into STL files by the CAD module, as shown in Figure 4a. Figure 4b shows the triangular facets in the STL file, which were described using three sets of X, Y and Z coordinates for each vertex and a normal unit vector indicating the direction of the facet. The three vertices were arranged counterclockwise. It must be mentioned that the triangular facets in the STL file were isolated and had no relation to the surrounding triangles. In other words, the STL file contained only graphical location information, not topological information. Therefore, in order to extract geometric information more effectively in the CAD model for subsequent path planning, it was necessary to pre-process STL files and create the relations among points, edges and facets [21]. In order to create the overall topological information at one time, the topology access should meet the following basic rules: Rule 1: Given a triangular facet, get its three edges and three vertices. Rule 2: Given an edge, get its two adjacent triangular facets and two vertices. Rule 3: Given a vertex, get all triangular facets in which the vertex is located. As shown in Figure 5, the respective data structures of vertices, edges and facets were designed according to the rules above. To store and read data more efficiently, three lists of vertices, edges and facets were created to store vertices, edges and facets respectively. The information of a vertex contained the index of the vertex, the coordinate of the vertex, the normal vector of the vertex and all indexes of the facet on which the vertex is located. The information of an edge contained the index of the edge, indexes of two end points of the edge and indexes of two facets on which the edge is located. The information of a facet contained the normal vector of the surface, the index of three vertices and the index of three edges. According to the characteristics of the format of an STL file, a topological relation reconstruction algorithm was designed to read and parse STL files. The information of each triangular facet was read one by one, including the vertex coordinates and normal vector. Each vertex, edge and facet was topologically related to the surrounding elements and stored in the corresponding lists, as in the following steps. As shown in Figure 5, the respective data structures of vertices, edges and facets were designed according to the rules above. To store and read data more efficiently, three lists of vertices, edges and facets were created to store vertices, edges and facets respectively. The information of a vertex contained the index of the vertex, the coordinate of the vertex, the normal vector of the vertex and all indexes of the facet on which the vertex is located. The information of an edge contained the index of the edge, indexes of two end points of the edge and indexes of two facets on which the edge is located. The information of a facet contained the normal vector of the surface, the index of three vertices and the index of three edges. As shown in Figure 5, the respective data structures of vertices, edges and facets were designed according to the rules above. To store and read data more efficiently, three lists of vertices, edges and facets were created to store vertices, edges and facets respectively. The information of a vertex contained the index of the vertex, the coordinate of the vertex, the normal vector of the vertex and all indexes of the facet on which the vertex is located. The information of an edge contained the index of the edge, indexes of two end points of the edge and indexes of two facets on which the edge is located. The information of a facet contained the normal vector of the surface, the index of three vertices and the index of three edges. According to the characteristics of the format of an STL file, a topological relation reconstruction algorithm was designed to read and parse STL files. The information of each triangular facet was read one by one, including the vertex coordinates and normal vector. Each vertex, edge and facet was topologically related to the surrounding elements and stored in the corresponding lists, as in the following steps. According to the characteristics of the format of an STL file, a topological relation reconstruction algorithm was designed to read and parse STL files. The information of each triangular facet was read one by one, including the vertex coordinates and normal vector. Each vertex, edge and facet was topologically related to the surrounding elements and stored in the corresponding lists, as in the following steps.
Step 1: Read in the three vertices (V 1 , V 2 , V 3 ) of a triangular facet in turn, and judge whether the vertex already exists in the list of vertices. If it does not exist, add it to the list of vertices; if it does, skip it.
Step 2: Read in the three edges (E 1 , E 2 , E 3 ) of the current triangular facet in turn, and judge whether the edge already exists in the list of edges. If it does not exist, the index of the two vertices that make up the edge will be written in this edge, and add this edge to the list of edges; if it does, skip it.
Step 3: When all vertices and edges of a triangular facet (F) are added to the corresponding lists, the index of the three vertices and three edges that make up the facet will be written in this facet. Then, add this facet to the list of facets.
Step 4: Repeat Step 1-3 until all triangular facets are read and recorded.
To illustrate this procedure further, the pseudocode of this procedure is also given as the following Algorithm 1. The local relations of a triangular mesh example extracted from the STL file are shown in Figure 6. The topological relation of edges after topological reconstruction and the quantity of the recorded data are listed in Tables 1 and 2, respectively. Comparing Figure 6a with Figure 6b, it was obvious that the data structure was significantly simplified after topological reconstruction. The data of vertices, edges and facets were recorded uniquely after topological reconstruction, which contained no redundant information and provided complete topological relation. The topological relations of the edges are listed in Table 2. According to such information, it was easy to find the two vertices of each edge and the two facets on which the edge is located. This advantage played a huge role in the subsequent path planning method.  In this section, the initial cladding path was generated from the CAD model based on the topological relations. A series of intersection point sets were obtained using a series of equidistant parallel planes intersecting the CAD model, which was called the slicing process. It was obvious that the process of intersection between the slicing planes and the CAD model was actually to find the intersection between the edges of the triangular facets and the cross planes. After all the intersection points were calculated, they were connected in a certain order to obtain the transversal lines.
In order to reduce the number of intersection calculations and avoid additional time consumption of sorting the intersection point sets, the slicing algorithm was designed based on topological relations of the edges. The core of the slicing algorithm was to constantly search the edge list for the edge that intersected with the current slicing plane and calculate the intersection point. The basic principle of the algorithm is, firstly, to traverse the list of edges to find an edge that intersects with the current slicing plane and, then, to search for next intersected edge in the adjacent   Table 2. Topology of edges after topological reconstruction.

The Slicing Algorithm Based on Topological Relation
In this section, the initial cladding path was generated from the CAD model based on the topological relations. A series of intersection point sets were obtained using a series of equidistant parallel planes intersecting the CAD model, which was called the slicing process. It was obvious that the process of intersection between the slicing planes and the CAD model was actually to find the intersection between the edges of the triangular facets and the cross planes. After all the intersection points were calculated, they were connected in a certain order to obtain the transversal lines.
In order to reduce the number of intersection calculations and avoid additional time consumption of sorting the intersection point sets, the slicing algorithm was designed based on topological relations of the edges. The core of the slicing algorithm was to constantly search the edge list for the edge that intersected with the current slicing plane and calculate the intersection point. The basic principle of the algorithm is, firstly, to traverse the list of edges to find an edge that intersects with the current slicing plane and, then, to search for next intersected edge in the adjacent edges of this edge according to the topological relation and to keep searching along the given direction until there are no intersected edges or return to the initial edge.
To distinguish the searched edge from the unsearched edge, a Boolean variable named isSearched was declared as an identifier in the data structure of the edge. The value of the identifier for all the edges was initially set as FALSE, which indicates none of the edges have been searched, and then changed to TRUE after being searched. The identifier was greatly useful for ensuring that each edge was accessed only once, which can significantly improve the efficiency of the slicing algorithm.
In addition, there were several important variables, namely the width, direction and range of the slicing process, which needed to be determined before the slicing algorithm was executed. The direction and width of the slicing process were specified by the user according to the specific laser cladding requirements. The range of the slicing process was determined by the maximum and minimum coordinates on the slicing direction in the list of vertices. When the slicing direction was determined, the slicing algorithm automatically calculated the range of the slicing process along the slicing direction.
Take the Z direction for example, the procedure of laser cladding path generation based on topological slicing is shown in Figure 7, where d is the slicing width, Z is the slicing direction and S is the value of the slicing plane in Z direction. The value of the slicing plane was initially set as the minimum value in Z direction and incremented by the value of the slicing width d. The slicing process started at the minimum coordinate and ended at the maximum coordinate.
Algorithms 2020, 13, x FOR PEER REVIEW 8 of 19 edges of this edge according to the topological relation and to keep searching along the given direction until there are no intersected edges or return to the initial edge.
To distinguish the searched edge from the unsearched edge, a Boolean variable named isSearched was declared as an identifier in the data structure of the edge. The value of the identifier for all the edges was initially set as FALSE, which indicates none of the edges have been searched, and then changed to TRUE after being searched. The identifier was greatly useful for ensuring that each edge was accessed only once, which can significantly improve the efficiency of the slicing algorithm.
In addition, there were several important variables, namely the width, direction and range of the slicing process, which needed to be determined before the slicing algorithm was executed. The direction and width of the slicing process were specified by the user according to the specific laser cladding requirements. The range of the slicing process was determined by the maximum and minimum coordinates on the slicing direction in the list of vertices. When the slicing direction was determined, the slicing algorithm automatically calculated the range of the slicing process along the slicing direction.
Take the Z direction for example, the procedure of laser cladding path generation based on topological slicing is shown in Figure 7, where d is the slicing width, Z is the slicing direction and S is the value of the slicing plane in Z direction. The value of the slicing plane was initially set as the minimum value in Z direction and incremented by the value of the slicing width d. The slicing process started at the minimum coordinate and ended at the maximum coordinate. Each single transversal line was generated in the following steps: Step 1: Find an edge that intersects with the current slicing plane but isSearched is FALSE. Figure out the intersection between this edge and the slicing plane and change isSearched of this edge to TRUE.
Step 2: Find all the adjacent edges of that edge according to the topological relationship between the points, edges and facets. If the isSearched of that edge is FALSE, figure out the intersection between that unsearched edge and the slicing plane; if it is TRUE, skip it.
Step 3: End the search when an unsearched edge cannot be found. Generate the path by joining all the intersections.
This slicing algorithm had two advantages. Firstly, through this method, it only needed to traverse the list of edges once to find the first edge intersecting the first current tangent plane. In the actual slicing process, only the edges contained in the triangular facet intersecting with the current tangent plane were searched and calculated, which greatly reduced the number of intersection Each single transversal line was generated in the following steps: Step 1: Find an edge that intersects with the current slicing plane but isSearched is FALSE. Figure out the intersection between this edge and the slicing plane and change isSearched of this edge to TRUE.
Step 2: Find all the adjacent edges of that edge according to the topological relationship between the points, edges and facets. If the isSearched of that edge is FALSE, figure out the intersection between that unsearched edge and the slicing plane; if it is TRUE, skip it.
Step 3: End the search when an unsearched edge cannot be found. Generate the path by joining all the intersections.
This slicing algorithm had two advantages. Firstly, through this method, it only needed to traverse the list of edges once to find the first edge intersecting the first current tangent plane. In the actual slicing process, only the edges contained in the triangular facet intersecting with the current tangent plane were searched and calculated, which greatly reduced the number of intersection calculations. Secondly, the topological slicing algorithm was automatically executed in a definite search direction. The ordered intersection sequence was directly obtained, which avoids the reordering and reconnection of the intersection set and greatly improves efficiency. In this way, the cladding path, which included a set of parallel, equidistant and directed transversal lines, was obtained directly.
To illustrate this procedure further, the pseudocode of this procedure is also given in Algorithm 2.

Simplification of the Intersection Data
It is difficult to directly interpolate curves when a curved surface part is processed using laser cladding. In engineering practice, the polyline, which includes multiple point-to-point line segments, is usually used to approximate the real curve line. The intersection points set directly obtained by the slicing algorithm was large and dense, which could not be directly used as the laser cladding processing point set. If it is directly used, the robot program file will be too large to execute, and the robot control system has to perform the interpolation point by point with extremely low efficiency. Furthermore, a greater degree of alloying and large dilution rate will occur during the laser cladding process when the cladding point set is too dense, which would pollute the composition of coatings and destroy its performance [22,23]. Therefore, on the premise of processing accuracy, it is essential to simplify the intersection point set and reasonably extract the effective processing points. The arch height error method [24] was introduced into this investigation as an effective mathematical method to optimize the discrete intersection point set.
The arch height calculation is shown in Figure 8. Known from the previous section, the set of discrete intersection points {P i } was already an ordered sequence. The allowable arch height error dε was the key parameter that affected the number and density of the final interpolation point set, which determined the accuracy and quality of the final cladding path. In this study, dε was usually set as the allowable maximum limit value of defocusing change according to specific laser cladding requirements.
Algorithms 2020, 13, x FOR PEER REVIEW 10 of 19 The arch height calculation is shown in Figure 8. Known from the previous section, the set of discrete intersection points {Pi} was already an ordered sequence. The allowable arch height error dε was the key parameter that affected the number and density of the final interpolation point set, which determined the accuracy and quality of the final cladding path. In this study, dε was usually set as the allowable maximum limit value of defocusing change according to specific laser cladding requirements. The general steps for finding interpolation points are described as the following: Step 1: Assuming that point Pi of the sequence is selected as the starting point, calculate the distance di+1 from the point Pi+1 to the line PiPi+2.
Step 2: If di+1 > dε, it means that the defocusing change exceeds the allowable range, then add Pi+1 as the valid interpolation point to the set of interpolation points. If di+1 < dε, it means that the defocusing change is within the allowable range, then go to the next step.
Step 3: Connect the Pi and Pi+3 points and calculate the distance from the points Pi+1 and Pi+2 to the line PiPi+3 as di+1 and di+2 respectively. The maximum dmax of the collection of {di+1, di+2} represents the maximum arch height of the arc segment PiPi+3. If dmax > dε, then add Pi+2 as the valid interpolation point to the set of interpolation points. If dmax < dε, continue to search forward until the maximum dmax of the collection of {di+1, di+2, di+3, ... di+n} exceeds the allowable arch height error dε, then add the point Pi+n as the valid interpolation point to the set.
Step 4: Repeat Step 1-3, while point Pi+n is selected as the new starting point. The simplification of the intersection point set usually started with the first point P1 in intersection points set {Pi} and ended with the last point Pm. Through the above method, all the effective processing interpolation points could be found and outputted as a sequence in the original order.

Calculation of the normal vector of cladding point
In order to facilitate the subsequent calculation of the position and pose of the robot TCP, the normal vector of vertex in the CAD model should be estimated first. The normal vector of each vertex was estimated by adding the normal vector of all triangles in which is located the center of gravity by area [25]. As shown in Figure 9, Aj denotes the area of the triangle, nj denotes the normal vector of each triangle, V denotes the common vertex, and N denotes the normal vector of the common vertex V. The normal vector N of vertex V can be estimated as follows: This formula indicated that the normal vector of the triangle facet with a smaller area contributed more to the normal vector N of the common vertex V. The general steps for finding interpolation points are described as the following: Step 1: Assuming that point P i of the sequence is selected as the starting point, calculate the distance d i +1 from the point P i+1 to the line P i P i+2 .
Step 2: If d i+1 > dε, it means that the defocusing change exceeds the allowable range, then add P i+1 as the valid interpolation point to the set of interpolation points. If d i+1 < dε, it means that the defocusing change is within the allowable range, then go to the next step.
Step 3: Connect the P i and P i+3 points and calculate the distance from the points P i+1 and P i+2 to the line P i P i+3 as d i+1 and d i+2 respectively. The maximum d max of the collection of {d i+1 , d i+2 } represents the maximum arch height of the arc segment P i P i+3 . If d max > dε, then add P i+2 as the valid interpolation point to the set of interpolation points. If d max < dε, continue to search forward until the maximum d max of the collection of {d i+1 , d i+2 , d i+3 , . . . d i+n } exceeds the allowable arch height error dε, then add the point P i+n as the valid interpolation point to the set.
Step 4: Repeat Step 1-3, while point P i+n is selected as the new starting point. The simplification of the intersection point set usually started with the first point P 1 in intersection points set {P i } and ended with the last point P m . Through the above method, all the effective processing interpolation points could be found and outputted as a sequence in the original order.

Calculation of the Normal Vector of Cladding Point
In order to facilitate the subsequent calculation of the position and pose of the robot TCP, the normal vector of vertex in the CAD model should be estimated first. The normal vector of each vertex was estimated by adding the normal vector of all triangles in which is located the center of gravity by area [25]. As shown in Figure 9, A j denotes the area of the triangle, n j denotes the normal vector of each triangle, V denotes the common vertex, and N denotes the normal vector of the common vertex V. The normal vector N of vertex V can be estimated as follows: This formula indicated that the normal vector of the triangle facet with a smaller area contributed more to the normal vector N of the common vertex V.
The normal vector of the processing interpolation point could be determined by the normal vectors of the two vertices on the edge, on which the interpolation point located. As shown in Figure 9b, P denotes an interpolation point, P 1 and P 2 denote two vertices on the edge containing P, and N, n 1 and n 2 denote the normal vectors of P, P 1 and P 2 , respectively. The normal vector N of P could be calculated using the following equation: where l 1 , l 2 and l denote the length of the line segment P 1 P, PP 2 and P 1 P 2 respectively. If the slicing plane crossed exactly one vertex, then the normal vector N of P was equal to the normal vector of the vertex.
Algorithms 2020, 13, x FOR PEER REVIEW 11 of 19 The normal vector of the processing interpolation point could be determined by the normal vectors of the two vertices on the edge, on which the interpolation point located. As shown in Figure  9b, P denotes an interpolation point, P1 and P2 denote two vertices on the edge containing P, and N, n1 and n2 denote the normal vectors of P, P1 and P2, respectively. The normal vector N of P could be calculated using the following equation: where l1, l2 and l denote the length of the line segment P1P, PP2 and P1P2 respectively. If the slicing plane crossed exactly one vertex, then the normal vector N of P was equal to the normal vector of the vertex.

Determination of the robot position and pose
In order to form a uniform and stable molten pool on the substrate during the laser cladding process, the laser beam is usually required to be perpendicular to the laser cladding surface. If the molten pool inclines during the cladding process, the shape and size of the coating section will be inconsistent, which affects the final quality of the coatings [26,27]. In this paper, therefore, the laser head was fixed on the worktable to keep it vertical all the time, and the workpiece clamped by the robot moved accordingly. In this mode, the laser beam was always perpendicular to the cladding surface, and the molten pool basically remained stable without tilting, which was beneficial for the formation of coatings. The where a = [ax ay az] T denotes the normal vector of the interpolation point, o = [ox oy oz] T is defined as the direction vector of the cladding processing path, and n = [nx ny nz] T is determined by n = o × a.
Matrix T represented the rotation matrix of the interpolation point relative to the workpiece coordinate system. However, in actual processing, the normal vector of the interpolation point was strictly required to coincide with the laser beam. Therefore, to meet this requirement, it was wise to set the interpolation point coordinate system {P} to be consistent with the world coordinate system {W} by adjusting the pose of the robot TCP. In other words, setting the pose of the robot made the vectors n, o and a coincide with the X, Y and Z axes of the world coordinate system, respectively, as shown in Figure 10.

Determination of the Robot Position and Pose
In order to form a uniform and stable molten pool on the substrate during the laser cladding process, the laser beam is usually required to be perpendicular to the laser cladding surface. If the molten pool inclines during the cladding process, the shape and size of the coating section will be inconsistent, which affects the final quality of the coatings [26,27]. In this paper, therefore, the laser head was fixed on the worktable to keep it vertical all the time, and the workpiece clamped by the robot moved accordingly. In this mode, the laser beam was always perpendicular to the cladding surface, and the molten pool basically remained stable without tilting, which was beneficial for the formation of coatings.
The coordinates and normal vector of the interpolation point were recorded as V = [Vx Vy Vz] T and a = [a x a y a z ] T respectively in the workpiece coordinate system. The position and pose of the TCP defined as P = [Px Py Pz α β γ] could be calculated using a series of transformations. The rotation matrix T composed of vectors n, o and a could be expressed as: where a = [a x a y a z ] T denotes the normal vector of the interpolation point, o = [o x o y o z ] T is defined as the direction vector of the cladding processing path, and n = [n x n y n z ] T is determined by n = o × a. Matrix T represented the rotation matrix of the interpolation point relative to the workpiece coordinate system. However, in actual processing, the normal vector of the interpolation point was strictly required to coincide with the laser beam. Therefore, to meet this requirement, it was wise to set the interpolation point coordinate system {P} to be consistent with the world coordinate system {W} by adjusting the pose of the robot TCP. In other words, setting the pose of the robot made the vectors n, o and a coincide with the X, Y and Z axes of the world coordinate system, respectively, as shown in Figure 10. In this case, the rotation matrix T' of the TCP coordinate system relative to the world coordinate system was equal to the rotation matrix of the workpiece coordinate system relative to the interpolation point coordinate system, expressed as the following: As shown in Figure 11, vector L represents the position vector of the laser focus in the world coordinate system; vector V represents the position vector of the interpolation point in the workpiece coordinate system; and vector D represents the laser defocusing quantity of cladding processing. The position vector V' of the interpolation point in the world coordinate system could be expressed as the following: ' ' Furthermore, the position coordinates of the robot TCP could be expressed as the following: In this case, the rotation matrix T of the TCP coordinate system relative to the world coordinate system was equal to the rotation matrix of the workpiece coordinate system relative to the interpolation point coordinate system, expressed as the following: As shown in Figure 11, vector L represents the position vector of the laser focus in the world coordinate system; vector V represents the position vector of the interpolation point in the workpiece coordinate system; and vector D represents the laser defocusing quantity of cladding processing. The position vector V' of the interpolation point in the world coordinate system could be expressed as the following: Furthermore, the position coordinates of the robot TCP could be expressed as the following: The pose coordinates of the robot TCP are described as RPY rotation angles: α (Roll angle), β (Pitch angle) and γ (Yaw angle) as shown in Figure 12, which respectively represented the rotation of the robot TCP coordinate system around the Z, Y and X axis of the world coordinate system.
By comparing Equation (4) and (7), a set of equations for α, β and γ could be described: By solving the equations, the expressions of α, β and γ could be obtained as follows: The pose coordinates of the robot TCP are described as RPY rotation angles: α (Roll angle), β (Pitch angle) and γ (Yaw angle) as shown in Figure 12, which respectively represented the rotation of the robot TCP coordinate system around the Z, Y and X axis of the world coordinate system. The pose coordinates of the robot TCP are described as RPY rotation angles: α (Roll angle), β (Pitch angle) and γ (Yaw angle) as shown in Figure 12, which respectively represented the rotation of the robot TCP coordinate system around the Z, Y and X axis of the world coordinate system.
By comparing Equation (4) and (7), a set of equations for α, β and γ could be described: By solving the equations, the expressions of α, β and γ could be obtained as follows: The compound rotation matrix RPY was: cos α cos β cos α sin β sin γ − sin α cos γ cos α sin β cos γ + sin α sin γ sin α cos β sin α sin β sin γ + cos α cos γ sin α sin β cos γ − cos α sin γ − sin β cos β sin γ cos β cos γ By comparing Equations (4) and (7), a set of equations for α, β and γ could be described: cos α cos β = n x sin α cos β = o x − sinβ = a x cos β sin γ = a y cos β cos γ = a z By solving the equations, the expressions of α, β and γ could be obtained as follows: According to the calculation results of Equations (6) and (9), the final position and pose of the robot TCP could be expressed as P = [Px Py Pz α β γ]. The robot positions and poses corresponding to all the interpolation points were obtained using the above method.

Automatic Path Generation Process
In order to realize the path planning method proposed in this paper, the authors developed a processing software based on C# and .net framework 4.5. The laser cladding path was generated and simulated using the processing software on a PC platform (CPU: Intel Core i3-4160 3.60 GHz; RAM 4.0 GB; OS: Microsoft Windows 10). An acetabular cup of an artificial hip joint with an outer diameter of 46 mm was selected as an example to verify the automatic generated path. The parameters related to path planning are listed in Table 3. The complete path generation process is shown in Figure 13. The STL model of the acetabular cup surface was sliced by five parallel planes in direction Z. Figure 13b shows five transversal lines with cladding spacing of 6 mm obtained by using the slicing algorithm based on a topological relation. Each single transversal line was connected by numerous intersection points in a laser scanning sequence. Figure 13c shows all the intersection points directly obtained. Figure 13d shows that the full path of the laser cladding consisted of five transversal lines end to end, each of which could be treated as single-pass cladding. To ensure the continuity of the cladding processing, the connecting lines were added between different single-pass cladding. When passing through the connecting lines, the robot moved without laser and powder. Figure 13e shows the direct result of the robot interpolation points with a defocus distance of 15 mm. Figure 13f shows the simplified result with an allowable arch height error of 0.5mm. Obviously, this not only guaranteed the shape precision of the surface machining, but also reduced the number of robot interpolations. Some key results of automatic path generation are listed in Table 4. Through topology reconstruction, the amount of data extracted and saved from the STL model was greatly reduced. The number of interpolation points were simplified from 850 to 74. The compression ratio of the interpolation points was approximately 8.7%, which contributed to reducing the size of the final robot program. On the other hand, the elapsed time of automatic path generation that involved software operations and computer processing was less than 30 seconds. When the full path was generated, the coordinates and normal vectors of the interpolation points were used to calculate the position and pose of the robot, which were finally transformed to the robot program. Some key results of automatic path generation are listed in Table 4. Through topology reconstruction, the amount of data extracted and saved from the STL model was greatly reduced. The number of interpolation points were simplified from 850 to 74. The compression ratio of the interpolation points was approximately 8.7%, which contributed to reducing the size of the final robot program. On the other hand, the elapsed time of automatic path generation that involved software operations and computer processing was less than 30 seconds. When the full path was generated, the coordinates and normal vectors of the interpolation points were used to calculate the position and pose of the robot, which were finally transformed to the robot program. In the same case, the manual teaching method needed more than two hours to edit the processing program. Therefore, compared with the traditional manual teaching method, the automatic path generation method proposed in this paper had great advantages in saving time, improving efficiency and ensuring safety.

Simulation Experiment
The robot model is NACHI MZ04 from Japan and the simulation module named "FD on Desk" matches with the robot. As shown in Figure 14, after adding the entry and exit points, the text file was imported into the simulation module and compiled into the executable file of the robot program. First of all, the simulation environment was set up according to the actual working conditions. Then, the compiled robot program file was inputted into the simulation model to run in automatic operation mode.  In the same case, the manual teaching method needed more than two hours to edit the processing program. Therefore, compared with the traditional manual teaching method, the automatic path generation method proposed in this paper had great advantages in saving time, improving efficiency and ensuring safety.

Simulation experiment
The robot model is NACHI MZ04 from Japan and the simulation module named "FD on Desk" matches with the robot. As shown in Figure 14, after adding the entry and exit points, the text file was imported into the simulation module and compiled into the executable file of the robot program. First of all, the simulation environment was set up according to the actual working conditions. Then, the compiled robot program file was inputted into the simulation model to run in automatic operation mode. Some different positions and poses of the robot during the simulation process are shown in Figure 15. After the automatic operation started, the TCP of the robot moved from the original position to the entry point and then experienced each interpolation point one by one along the path. Finally, it went through the exit point and returned to the original position. The results showed that the TCP of the robot moved exactly along the planned interpolation path. Furthermore, in the whole simulation process, the robot achieved a complete and continuous movement without mechanical collision occurring or exceeding the range of motion. The relative position relationship between the laser beam and the surface of the acetabular cup was shown through different angles. During the laser cladding process, the laser beam was always vertical and perpendicular to the processed surface. Some different positions and poses of the robot during the simulation process are shown in Figure 15. After the automatic operation started, the TCP of the robot moved from the original position to the entry point and then experienced each interpolation point one by one along the path. Finally, it went through the exit point and returned to the original position. The results showed that the TCP of the robot moved exactly along the planned interpolation path. Furthermore, in the whole simulation process, the robot achieved a complete and continuous movement without mechanical collision occurring or exceeding the range of motion. The relative position relationship between the laser beam and the surface of the acetabular cup was shown through different angles. During the laser cladding process, the laser beam was always vertical and perpendicular to the processed surface.  Figure 16 shows the changes of the robot pose angles during the simulation process. From those pictures, we can see the roll angle α changes from −180° to 180°, which indicates the robot TCP rotates around the X axis of the world coordinate system for circular motion, and the pitch angle β and the yaw angle γ changes from −90° to 90°, which indicates the robot TCP rotates around the Y and Z axis of world coordinate system in a half circular motion, respectively. The data results are consistent with the desired results of robot motion. The five cycles of smooth robot pose angles indicates that the robot TCP completes the whole motion correctly along the full path consisting of five parallel transversal lines. There are no singularities and drastic change of the robot pose in the motion process. In a word, the simulation results preliminarily confirmed the correctness of the path automatically generated from the path planning module as well as the path planning method proposed in this paper.

Conclusions
In this paper, the laser cladding path planning for artificial joints was mainly investigated. In order to improve the efficiency of laser cladding on the surface of artificial joints, a path planning method for a laser cladding robot on an artificial joint surface based on topology reconstruction and a novel processing mode were designed. The main conclusions are as follows:  Figure 16 shows the changes of the robot pose angles during the simulation process. From those pictures, we can see the roll angle α changes from −180 • to 180 • , which indicates the robot TCP rotates around the X axis of the world coordinate system for circular motion, and the pitch angle β and the yaw angle γ changes from −90 • to 90 • , which indicates the robot TCP rotates around the Y and Z axis of world coordinate system in a half circular motion, respectively. The data results are consistent with the desired results of robot motion. The five cycles of smooth robot pose angles indicates that the robot TCP completes the whole motion correctly along the full path consisting of five parallel transversal lines. There are no singularities and drastic change of the robot pose in the motion process.  Figure 16 shows the changes of the robot pose angles during the simulation process. From those pictures, we can see the roll angle α changes from −180° to 180°, which indicates the robot TCP rotates around the X axis of the world coordinate system for circular motion, and the pitch angle β and the yaw angle γ changes from −90° to 90°, which indicates the robot TCP rotates around the Y and Z axis of world coordinate system in a half circular motion, respectively. The data results are consistent with the desired results of robot motion. The five cycles of smooth robot pose angles indicates that the robot TCP completes the whole motion correctly along the full path consisting of five parallel transversal lines. There are no singularities and drastic change of the robot pose in the motion process. In a word, the simulation results preliminarily confirmed the correctness of the path automatically generated from the path planning module as well as the path planning method proposed in this paper.

Conclusions
In this paper, the laser cladding path planning for artificial joints was mainly investigated. In order to improve the efficiency of laser cladding on the surface of artificial joints, a path planning method for a laser cladding robot on an artificial joint surface based on topology reconstruction and a novel processing mode were designed. The main conclusions are as follows: In a word, the simulation results preliminarily confirmed the correctness of the path automatically generated from the path planning module as well as the path planning method proposed in this paper.

Conclusions
In this paper, the laser cladding path planning for artificial joints was mainly investigated. In order to improve the efficiency of laser cladding on the surface of artificial joints, a path planning method for a laser cladding robot on an artificial joint surface based on topology reconstruction and a novel processing mode were designed. The main conclusions are as follows: 1. Different from traditional laser cladding mode, the presented mode of the laser beam fixation during workpiece movement met the requirements of the laser cladding process of the artificial joints and ensured the safety of the operator and equipment.
2. The path planning method based on the topological relations greatly reduced the amount of intersection calculation and directly obtained the sequential and simplified cladding path, which can be translated to the robot program with a small scale.
3. The planned path on the surface of a 46 mm artificial acetabular cup model was carried out by the simulation experiment. The results show that the robot accurately drove the artificial acetabular cup to move along the planned cladding path without mechanical collisions, trajectory deviations and so on.
Many facets of path planning remain for future investigation. We will still focus on the efficiency and quality of the path planning method. Moreover, we will combine laser cladding technology to investigate the quality of the coatings processed by applying the path planning method.