A fast methodology for generating skeletal FEM with detailed human geometric features based on CPD and RBF algorithms

Due to the significant effects of the human anatomical characteristics on the injury mechanism of passenger in traffic accidents, it is necessary to develop human body FEM (Finite Element Model) with detailed anatomical characteristics. However, traditional development of a human body FEM is an extremely complicated process. In particular, the meshing of human body is a huge and time-consuming project. In this paper, a new fast methodology based on CPD (Coherent Point Drift) and RBF (Radial Basis Function) was proposed to achieve the rapid developing the FEM of human bone with detailed anatomical characteristics. In this methodology, the mesh morphing technology based the RBF was used to generate FEM mesh in the geometry extracted from the target CT (Computed Tomography) data. In order to further improve the accuracy and speed of mesh morphing, the target geometric feature points required in the mesh morphing process were realized via the rapid and automatic generation based on the point-cloud registration technology of the CPD algorithm. Finally, this new methodology was used to generate a 3-year-old ribcage FEM consisting of a total of 27,728 elements with mesh size 3–5 mm based on the THUMS (Total Human Model for Safety) adult model. In the entire process of generating this new ribcage model, it only took about 2.7 s. The average error between the new FEM and target geometries was only about 2.7 mm. This indicated that the new FEM well described the detailed anatomical characteristics of target geometry, thus importantly revealing that the mesh quality of the new FEM was basically similar to that of source FEM.

Passenger injury caused by traffic accidents is a serious public health issue worldwide 1,2 . Injury mechanism of passenger resulting from different traffic accidents plays an important role in the theoretical basis for solving such public health issue. It has been demonstrated that the human anatomical characteristics have significantly effects on the injury mechanism of passenger 3 . For example, Ridella et al. 4 reported that obese elderly and child passengers were more likely to be injured than those with normal body characteristics. Therefore, it is greatly significant to study the influences of human anatomical characteristics on injury mechanisms towards the protection of special passengers.
As of today, computational simulation has become one of the main methods to study injury mechanism and establish injury tolerances 5 . In particular, the detailed biomechanical responses of human tissue pertaining to injury severity and location, such as strain and stress, can be predicted by the human body FEM (finite element model). In addition, the human body FEM can accurately characterize the anatomical features of human body. The human body FEM has become one of the most widely used human injury assessment tool in the field of vehicle safety 6  www.nature.com/scientificreports/ usually includes: 1. creation of human body geometric model from CT (Computed tomography) and MRI (Magnetic resonance imaging) data, 2. FE meshing of human body geometric model, and 3. boundaries, loading, and verification of model. One of the most time-consuming is the meshing process of the geometric model. Existing meshing methods mainly include Delaunay, Advancing Front Technique, Mapping, Sweeping 11,12 . In practical use, the geometric model is decomposed and the meshes are generated mainly through manual interaction. However, for the human body structure with very complex geometric details, the mesh generated by the above methods is of poor quality and can hardly meet the analysis requirements. It requires experienced researchers to improve mesh quality by Laplacian smoothing, elements topology optimization and other operations 13,14 .
Due to the complexity of human body FEM development, the method to obtain new FEM based on the existing basic model through mesh deformation technology has been widely developed. The scaling method was firstly proposed. For example, Vavalle et al. 15 and Schoell et al. 16 scaled the 50th male FEM in the GHBMC to obtain a 95th adult male and a 65-year-old male FEM, respectively. In the scaling method, the new FEM is usually obtained by scaling the body parts of the existing basic FEM with different ratios without need of the detailed geometric data of the target geometry. This advantage makes the scaling method widely used in literatures [17][18][19] . However, the disadvantage of the scaling method is also obvious in aspect that the detailed geometric differences between the body parts of the target FEM and the existing basic FEM are not reflected. Considering this disadvantage of the scaling method, the UMTRI (University of Michigan Transportation Research Institute) and Hunan University recently proposed a mesh morphing method based on feature points and RBF (Radial Basis Functions) [20][21][22][23] . Mesh morphing is the smooth transition of a FEM into another similar FEM, where the first model is called the source FEM and the second is called the target FEM. Mesh morphing usually consists of three steps: Firstly, a large number of corresponding feature points are selected at appropriate locations of source FEM and target geometry. Then, the corresponding relationship between them is established through feature points, such as the corresponding relationship between vertices, edges and surfaces. Finally, the meshes of source FEM are mapped to the target geometry by RBF, so the target FEM is obtained. The target FEM generated by mesh morphing can retain the detailed geometric features of target geometry well. However, the limitation of mesh morphing is that a large number of feature points need to be selected manually, which is very time-consuming and laborious. For example, using mesh morphing to generate a ribcage FEM with 27,728 elements usually needs to manually select more than 1,000 feature points 24 . Moreover, once the sequence and number of feature points on the source FEM and the target geometry are inconsistent, the process of mesh morphing cannot be carried out. Considering all of these is necessary to improve such a time-consuming, labor-intensive, and error-prone step.
In this paper, in order to avoid the disadvantage of the mesh morphing method, an automatic generating feature points method using the CPD (Coherent Point Drift) algorithm was proposed. This new methodology was then applied to automatically generate feature points for different human bones, such as ribcage, pelvis, humerus, radius, tibia, and ulna. Furthermore, using the generated feature points, the FE meshing of these human bones was generated by the mesh morphing with RBF. In these applications, generating about 200 and 2000 feature points only takes about 2 s and 24 s, respectively. The quality of FE mesh obtained by using the automatically generating feature points is basically the same as that of FE mesh prior to the morphing. Results reveal that this method capable of generating the feature points automatically is faster and more accurate than the manual extraction.

Method
The Ethical Committee of Shanghai Ninth People's Hospital approved this retrospective study. And written informed consent was obtained from all the participants. All methods were performed in accordance with relevant named guidelines and regulations.
A whole process from CT data to the FEM using the fast mesh morphing method is shown in Fig. 1. First of all, the geometry of target new model described by point-cloud was extracted from CT data (ribcage, pelvis, humerus, radius, tibia and ulna). In this process, the ribcage was taken as an example of the target new model. Second, the corresponding source FEM of ribcage was split from THUMS, and outermost mesh nodes were extracted as the source point-cloud. Third, a rough registration between the target and source point-clouds was conducted through PCA (Principal Component Analysis) with CPD algorithm to obtain feature points. Finally, the feature points were used to morph the source FEM for target FEM.
The geometry of target new model. In this study, a CT data of the whole-body for a 3-year-old male child was used to develop the target shape. By adjusting the CT value, the data related to the bone was extracted. Then, the data of bone was thresholding, and smoothed. Sometimes, the CT image quality is not good enough to capture all the geometrical details. For this, threshold adjustment was used to improve the completeness of the geometry. The target shape was then repaired by editing mask. Generally speaking, the complete usable shape mask can be extracted through the preprocessing process described above. The pre-processing process is shown in Fig. 2. The data processed above were exported in the form of a point-cloud.
The source FEM. The skeletal FEMs were split from THUMS and marked as FEM S for source FEMs. Then the outermost nodes of each FEM S were extracted as the source point-cloud for each bone and marked as P S .
Automatically generating feature points. Figure 3 shows the method overview for automatically generating feature points. Firstly, the target point-cloud extracted from CT data was filtered. Then the coordinate system of target point-cloud and source point-cloud was unified through rough registration. Finally, non-rigid registration was carried out from source point-cloud to target point-cloud to obtain target feature points based on CPD algorithm. www.nature.com/scientificreports/ Radius filtering of point-clouds. When the point-cloud of target shape was generated from CT data, it usually had a large number of points with outliers (i.e., noise points) in Fig. 3a. These noise points may cause the mismatches during the registration of the target and source point-clouds. Therefore, the radius filter was adopted to remove these noise points. In the area with a point as the center and d as the radius, if the number of points is less than K, the center point will be removed by the radius filter. In this study, the radius d and the number of points K was defined as 2 mm and 5, respectively. As shown in Fig. 3b, the noise points in the point-cloud of target shape were significantly reduced by the radius filter. The point-cloud of target shape without noise points was marked as P C_I .   www.nature.com/scientificreports/ Coordinate system unification. Usually, there is a huge difference of spatial position between P C_I and P S , because these two point-clouds were obtained in different coordinate systems. Therefore, it is necessary to unify the coordinate systems of these two point-clouds through matrix operations in Eq. (1).
where P C is the transformed point-cloud of P C_I in the coordinate system of P S ; R 0 is rotation matrix; T 0 is translation matrix. Using Eq. (1), the difficult step is to find out the rotation and translation matrices. Considering this, the PCA-based coarse registration method was adopted in this study, as shown in Fig. 3c 25 . In the PCA-based coarse registration, principal component analysis was used to reveal the main distribution direction of point-cloud and reduce dimension of the data. Therefore, the PCA-based coarse registration method is mainly based on the global principal axis direction of point-cloud data for registration. Firstly, the covariance matrix of P C_I and P S was calculated. And then, the main feature component, namely the global principal axis direction of the pointcloud data, was calculated according to the covariance matrix to make up the principal axis direction matrices. So, R 0 can be obtained by the principal axis direction matrices by Eq. (2). Finally, the coordinate values of center point calculated from two point-clouds and rotation matrix R 0 were used to obtain the translation matrix T 0 using Eq. (3).
where U X and U P are the 3*3 principal axis direction matrices of point-clouds P S and P C_I , respectively; P S and P C_I are the coordinate values of center point of point-clouds P S and P C_I , respectively.
Non-rigid registration. The P C and P S were aligned on the principal axis directions after the Coordinate System Unification. However, the mesh morphing technology is based on the same number of source and target feature points to achieve element mapping transformation. Hence it is necessary to create an equal number of source and target feature points in P C and P S . Considering the P S as extracted from the outermost nodes of the FEM with high quality mesh, therefore, the P S can be directly used as the source feature points. For the target feature points, a non-rigid registration was adapted to mapping the source feature points to the P C . The target feature points generated by this method not only can be consistent with the source feature points in number, but also have a similar distribution position in human body geometry to the source feature points. This can effectively reduce the distortion of the generating element in the mesh morphing. This non-rigid registration named as alignment depicted in Fig. 3d was realized by the CPD algorithm 26 .
The essence of obtaining target feature points by the non-rigid registration based on the CPD algorithm is to find out an accurate transformation matrix marked as T . Using this accurate T to transform P S can obtain a new point set P c_c which deems to be as similar as possible to point-cloud P C . P c_c can be used as the target feature points in the mesh morphing technology. In order to obtain an accurate T , the Gaussian Mixture Model (GMM) was adapted to address this problem in CPD algorithm. In the GMM, P S was considered as the GMM centroids and P C was considered as the GMM generated point-cloud. In other words, P S was regarded as a correct standard point-cloud, and P C was the point-cloud composed of many scattered points around P S . The relationship set between P S and P C in the GMM can be expressed by Eqs. (4) and (5). From Eq. (4), it can be found that the probability of the existence of a point in P C was described as the sum of the distance between this point and each GMM centroid (each point in P S ). Since P c_c obtained by transformation of P S through the accurate T was as completely coincident with P C as possible. It should also be pointed out that there is always a point in P c_c that can coincide with the corresponding point in P C (the distance between these two points is zero). Therefore, different T can be tried repeatedly to transform P S to obtain different P c_c . The possibility of all point in each P c_c was summed up by Eq. (6). This cumulative sum obtained by Eq. (6) can be used as the rating of the different T : the greater of the cumulative sum implies the higher score and more accurate of T . Accordingly, finding the accurate T needs to calculate the maximum value of the function described by Eq. (6). By maximizing the likelihood function described in Eq. (6), the parameters ( θ and σ 2 ) in T can be obtained. Finally, P c_c were calculated from transformation of P S by T.
where p(P C ) is represented as the generation probability of each point in P C ; N is the number of points in P C ; M is the number of points in P S , and w is the weight factor; p( P C |m) is the probability of the points in P C generated by each point in P S and can be calculated by Eq. (5).
(1) www.nature.com/scientificreports/ where θ is the parameter set consisting of the rotation matrix R , translation matrix T , and deformation matrix X ; Y m is the center of GMM model; X is the points generated by the GMM model; D is the dimensions of points in P C and P S .

Mesh morphing technology.
After the alignment of P C and P S in Section "Non-rigid registration", the outermost nodes (point-cloud P S ) in the source FEM were already mapped to the target shape described by P C . If the internal nodes in the source FEM can be mapped regularly into the target shape as expressed in Eq. (7), implying that target FEM with detailed target geometric characteristics will be obtained. This method for obtaining a new FEM is called as the mesh morphing technology. In 2016, Wang et al. 24 firstly proposed a mesh morphing technology based on RBF with kernel function TPS (Thin Plate Spline) 27 expressed as given in Eq. (8). In the RBF with kernel function TPS as given in Eq. (9), the source and target feature points obtained in Section "Automatically generating feature points" was used as the control points to calculate the weight coefficients in Eq. (8) by Eq. (10), then the internal element nodes associated with control points are smoothly transformed from the source FEM to the target FEM by Eq. (7).
where f x, y, z is the transformation function for mapping the internal nodes in the source FEM into the geometry of target new model; x where p x, y, z is a low order polynomial; ∅ is the kernel function representing the TPS in this paper; w i is the weight coefficient, is the coordinate values of i node in source feature points; h is representing the number of source feature points.
where W = [w 1 , w 2 , w 3 , . . . , w h ] ; Q is represented the coordinate values of the source feature points; V is represented as the coordinate values of the target feature points; A = a 0 , a x , a y , a z ; K is the distance between each source and its corresponding target feature point.

Results
In order to verify the method proposed in Section "Method", six skeletal FEMs separated from the adult 50 th THUMS including ribcage, pelvis, humerus, radius, tibia, and ulna were successfully morphed into the corresponding target shape extracted from a 3-year-old male child CT data as shown in Table 1 Geometric error between target shape and FEM. The distance between the outermost mesh nodes of target FEM and corresponding points of target point-cloud was used to evaluate the geometric error between target shape and FEM as shown in Fig. 4. It can be found that the average geometric error of each model is less than 3 mm, and particularly the average geometric errors of the humerus, radius, tibia and ulna model are less than 1.5 mm. Even the maximum geometric errors of the humerus, radius, tibia and ulna model are less than 5 mm. However, the maximum geometric errors of ribcage and pelvis FEMs are 15.332 and 14.645 mm, respectively. This is mainly due to the filtering of ribcage and pelvis original target point-clouds. The geometric features of ribcage and pelvis are complex, especially their original target point-clouds extracted from CT exhibits a large number of detailed geometric features, thus easily leading to incorrect point-cloud registration. Therefore, their original point-clouds are filtered to reduce useless detailed geometric features. However, the geometric error is the comparison between the target FEM and the original target point-cloud, showing some parts with detailed geometric features are relatively large.
Mesh quality. The mesh quality of solid elements in target FEMs including Jacobian, Warpage, Skew, Aspect ratio et al. were checked and listed in Table 2, showing that the overall mesh quality of the target FEMs was good and basically meets the requirements of finite element analysis. It was generally considered acceptable if the minimum Jacobian of the mesh is ≥ 0.2. However, it should be noted that the mesh quality of solid elements in target FEMs was a slightly inferior than that of solid elements in source FEMs. For example, compared to the source FEMs, the minimum Jacobians of the target humerus and ulna models were decreased from 0.39 and 0.31 to 0.25 and 0.22, respectively. Because the target feature points generated in Section "Mesh morphing technology" described the detailed geometric features of the target shape as much as possible, this led to a less uniformity of the elements distribution with these target feature points as the outermost nodes and the poor element quality. The mesh quality degraded by this factor can be improved with mesh smoothing.

Discussion
In the process of mesh morphing, it was found that the target feature points have significant influences on target FEMs. In this section, the influences of the generation methods and number of target feature points on the geometric error and mesh quality of target FEMs will be discussed.
The influences of generation method of feature points. In this study, the non-rigid registration based on the CPD algorithm was adapted in Section "Automatically generating feature points" to obtain the target feature points. In the non-rigid registration, using coordinate translation, rotation, scaling and local deformation to generate target feature points from source feature points can better describe the geometric characteristics of the target point-cloud. In fact, the target feature points can also be obtained by rigid registration. For CPD algorithm, rigid registration is a transformation that does not change the relative position between two points in the point-cloud, including translation, rotation and scaling. Non-rigid registration can change the relative position between two points in the point-cloud. It can cause local deformation by nonlinear matrix. We did both registrations by changing the θ in Eq. (6). The θ in rigid registration include rotation matrix R , translation matrix T and scaling matrix S . Instead of rigid registration, the θ contains an additional nonlinear matrix N. As shown in Fig. 5, when using non-rigid registration, ribcage, pelvis, humerus, radius, tibia, and ulna FEMs generation required 44.61 s, 21.9 s, 8.192 s, 5.581 s, 5.988 s, and 4.933 s respectively. However, due to the rigid    Fig. 6, the geometric errors including maximum and average errors of target FEMs generated by the nonrigid registration method are lower than that generated by the rigid registration. This is mainly because the nonlinear transformation is not adapted to make the target feature points closer to the geometric features of the target shape in the rigid registration method. Especially for humerus, radius, tibia, and ulna, their detailed geometric characteristics of source and target point-clouds differ greatly. Therefore, compared with target FEMs generated by the rigid registration, the average geometric error of target FEMs generated by the non-rigid registration are decreased by 73.3%, 72.2%, 77.7% and 66.5%, respectively. Figure 7 shows that the mesh quality of target FEMs generated by rigid registration are similar to that of source FEMs shown in Table 2. This is mainly because the translation, rotation and scaling in the rigid registration method did not cause local deformation of the mesh elements to affect the mesh quality of target FEMs. Therefore, the mesh quality of target FEMs generated by rigid registration were also higher than those of the corresponding target FEMs generated by non-rigid registration.
The influences of number of target feature points. Different number of target feature points have different functional capability to describe the geometric features of target point-cloud, affecting the generation of nodes for the target FEMs in the mesh morphing process. Therefore, the influences of number of target feature points were described in this section. Selecting 100%, 50%, and 10% of source features points by the uniform sampling method were respectively used to generate the corresponding number of target feature points by nonrigid registration for further generating target FEMs. It can be seen from Fig. 5, as the number of feature points decreases, the generation time of target FEMs decreases gradually. When 10% of source features points are selected, the ribcage FEM with 27,728 elements needs only 2.7 s. Figure 8 shows the geometric errors including the maximum and average errors of target FEMs increased with decreasing of the number of target feature points. The influences of number of target feature points on the geometric errors are slightly in the target FEMs of ribcage and pelvis, but significantly in the target FEMs of humerus, radius, tibia, and ulna. This is primarily due to the detailed geometric characteristics of their source and target point-clouds differ significantly. Therefore, more target feature points are necessary to describe the detailed geometric characteristics of target point-clouds.
Data shown in Table 3, reveal the mesh quality of the target FEMs generated by the small number of target feature points is better. This is mainly due to fewer target feature points that ignore some detailed geometric features of the garget point-cloud to reduce the local deformation of mesh elements. Similarly, in the general mesh generation, it is also important and necessary to balance the mesh quality and the details of geometric features.

Limitations.
Through presentation and discussion of the method proposed in this study, the source FEMs can be quickly morphed to the target FEMs more consistent with the geometric features of the target pointcloud. However, in order to describe more detailed geometric feature of target point-cloud, a non-rigid registration method is used to locally deform the source feature points to generate target feature points. This local deformation caused a decrease in the mesh quality of the target FEM compared to that of the corresponding source FEM. Especially when the detailed geometric features of the source and target point-clouds differ significantly, the mesh quality of the generated target FEMs decreases more significantly. For the meshes with acceptable quality, we usually improve them by smoothing, elements optimization, nodes optimization and so on. Therefore, it is necessary to conduct more in-depth study on how to balance the mesh quality and the details of geometric

Conclusion
In this study, a fast-morphing methodology for generating human skeletal FEM with detailed geometric features based on CPD and RBF algorithms was proposed. This is indeed to be the first time that the CPD algorithm is applied to the traditional mesh morphing technology. Using this algorithm enables realization of the www.nature.com/scientificreports/ improvement from manually placing feature points to automatically generating feature points. In addition, this method can directly use the point-cloud of the target geometric data for finite element modeling to save a lot of work on the reverse modeling. This fast-morphing methodology was successfully used to morph several human skeletal FEMs extracted from THUMS adult model to 3-year-old child skeletal FEMs. The morphing results proved that the human skeletal FEMs generated by this fast-morphing methodology has resulted in small geometric errors and high mesh quality. Limitations on this approach will be continued to be investigated further in research pertaining to this area in the future. www.nature.com/scientificreports/

Data availability
The datasets used and analyzed during the current study available from the corresponding author on reasonable request.