A Corresponding Method of 3D Skull Feature Points Based on Voxel Models and Multi-Geometric Feature Constraint

In the craniofacial reconstruction method based on statistical model, several feature points need to be calibrated on each skull. Due to the large number of feature points, the order of calibration is difficult to orderly, so it is necessary to carry out point correspondence to the calibrated feature points, and it is convenient to carry out subsequent statistical analysis. Therefore, a 3D skull correspondence algorithm based on voxel models and multi-geometric feature constraint is proposed in this paper. First, the skull model is unified into the Frankfurt coordinate system. According to the corresponding relationship between the reference skull and the corresponding set of skull points, the matching skull is deformed by TPS, and the skull is coarsely registered; Then, a voxel model of reference skull and skull to be matched is constructed to solve multi-geometric features of feature points set of the model; Finally, under the joint constraint of multi-geometric features, the best counterpart is achieved. Experimental results show that the proposed algorithm improves the accuracy and efficiency significantly.


Introduction
In the craniofacial restoration based on statistical models, the point correspondence of craniofacial data is an indispensable step in the establishment of craniofacial statistical model, and the accuracy of the corresponding results directly affects the effect of craniofacial restoration. The skull point correspondence refers to the mapping relation between the points of different skull models with the same physiological structure. As the skull model complex structure in a point corresponding relationship, the need for the following requirements of the skull model: skull model with equal number of vertices; vertex position the same number of different models with physiological consistency, such as tenth point sample A is the point between the eyebrows, is a sample of tenth point B is the point between the eyebrows the same location; samples need registration to space, avoid the influence of affine transformation.
At present, scholars at home and abroad have little research on the specific problem of 3D skull feature correspondence. The first point is the method of manual marking corresponding point, Cootes et al. [1] defines the key point of the structure, to establish the corresponding relationship through the same feature points manually mark position, the number of 3D skull model is very much, the manual alignment is very time-consuming; Chen et al. [2] proposed the nearest point correspondence method based on Euclidean distance. However, when the relative position gap between the two models is large, the corresponding results are unstable and the time complexity of the algorithm is increased; Marani et al. [3] introduced the concept of deleting masks and improved the iterative nearest point (ICP) algorithm. Although the method can reduce the error, the hypothesis is difficult to meet or satisfy in many times, and the cost of iterative computation is very large in the search process of the nearest point; Duan et al. [4,5] ICP through a coarse registration random sampling points as feature points corresponding to the TPS, then the TPS deformation, and the iteration times; Zhang et al. [6] are the corresponding points of anatomy after ICP coarse registration, which are the characteristic points of TPS. Because TPS is a global deformation algorithm, every TPS distortion will affect the result of the previous registration; Deng et al. [7] uses a radial basis function with tight support instead of TPS for accurate registration in ICP rough registration, which only adjusts the local area with large registration error; Du et al. [8] proposed the 3D skull point correspondence method based on ICP and TPS. The method introduced the spherical volume integral in-variants of multi-scale constraints to realize the corresponding point set of the reference skull and the point of the target skull. Liu et al. [9] proposed an improved ICP 3D point cloud registration method. This method can improve the speed and accuracy of point cloud registration and the stability. However, the large point cloud data need to be further studied and improved in the search strategy of the corresponding point. Therefore, a 3D skull correspondence algorithm based on voxel models multi-geometric feature constraint is proposed in this paper. The main steps of the method are: 1)The reference skull and the corresponding skull are unified into the Frankfurt coordinate system, and the corresponding skull model is deformed by TPS, so that the shape of the corresponding skull is approximately coincident with the reference skull shape; 2)Includes the establishment of vertices with one order adjacency voxel model for skull model, the calculation of differential properties reference vertices in skull model: the average normal and Gauss curvature, and according to the size of the Vertex Curvature do in descending order；3)In order to find the best matching points between the corresponding skull and the reference skull, the vertex of the curvature in the skull is selected as the corresponding point. The algorithm flow diagram is shown in Figure 1:

Matching skull
Reference skull Frankfurt coordinate system normalization; The corresponding skull TPS is deformed, and the coarse registration is completed To establish the voxel model of the skull; The normal vector and Gauss curvature of the vertex are calculated According to the order of curvature, the corresponding vertices are selected in turn, and the searching range of the vertex in the corresponding skull is determined According to the multi geometric features, the corresponding correspondence of the best points is established

Skull Rough Registration Based on TPS
In this paper, the method of literature [10] is used to match the skull model, because a linear polynomial and a thin plate spline are usually attached to the radial basis function when the thin slice is applied to the radial basis function.
Among them, the linear polynomial of affine transformation is ( ) x f , the thin plate spline is 2 (r) r lg(r) ϕ = , i r x x = − , and the corresponding linear function of the basis function is exactly corresponding to the affine component of the transformation between the feature points. Therefore, the slice spline function can describe the shape feature of feature points distribution, and has affine invariance. In this paper, the method can be used to establish a smooth interpolation function of the surface of the skull model, which can realize the non-rigid deformation and approximately coincide between the models. As shown in Figure 2(a), a matching result between the skull and the skull to be matched to the same coordinate system is shown. Figure 2

To Establish the Voxel Model of the Skull
The voxelization refers to the geometric representation of 3D skull model into the model closest to the voxel representation, form data set, which contains not only the surface information of the model, but also describes the internal attribute of the model, similar to the two-dimensional pixel, extended from the two dimension to the cube unit three.
The specific steps of the voxel method for skull model implementation are as follows: (1) According to the vertex data of the reading, the maximum value of max X , max Y , max Z and minimum min X , min Y , min Z on the three axes of X , Y and Z in the reference skull are calculated respectively.
(2) According to the maximum value obtained by the first step and the minimum value is kρ , the Here refers to the element length of dough mesh unit, according to the experimental results, when at 12-27, corresponding to the point of higher accuracy, each voxel usually contains only one or no vertex, the utility model has the advantages of small volume of voxels, including the vertex space description ability. From this we can see that voxel is a cube of size Among them, x, y, z domain is an array, save the coordinates of the current voxel contains vertices (may be more than one point), voxNum domain is an array, save the voxels contained in the vertex index, adjVer domain for cellular array, each line is only one cell line number, save the corresponding voxNum array a neighbor, initialized to the empty; (4) The coordinates of each vertex of the template was X , Y and Z three direction unit number, such as the current coordinates (x, y, z), the ordinal number of voxel in the X direction ). Index the current point, serial number, and its coordinates into voxel units (i,j,k); (5) According to the read skull data, all the two vertex numbers including the current vertex (x, y, z) are stored in the adjVer domain of the voxel index , and the first neighborhood relation of each vertex is established.
Through the above steps structure of the voxel model has the advantages of simple operation, can order neighborhood information convenient access to the vertices can be effectively constrained to the corresponding vertex in the reference skull matching search scope of the corresponding point in the skull, reducing the time complexity of the algorithm, the most important is the removal of the skull reference point corresponding process in the current vertex in the matching point in the skull not avoided, serious mismatches.

Multi-geometric Feature Extraction of Skull
Because the geometry of skull model is complex and diverse, the traditional surface registration and point correspondence method can not satisfy the accurate establishment of craniofacial correspondence. For the complex geometric model, the vertex differential attribute often has a strong representation capability, so we introduced vertex differential properties of the skull model, to determine the relation between the vertices with differential attribute weighted distance and Euclidean distance constraints between vertices, overcome the defect of pure Euclidean distance decision point corresponding relationship is not accurate enough.

The calculation of the normal vector of vertices
The calculation model of the skull vector. There are many methods, according to the first order neighborhood information obtained during the voxelization of vertices, the weighted average method on the adjacent triangle area and the angle of the vector method with first order adjacency domain current vertex i v to calculate the normal vector, the first-order neighborhood relationship using voxel model about the current vertex, to calculate the effective constraint method. The first order neighborhood of vertex i v is shown in Figure 3. The mean vector of the vertex i v is calculated as follows:

Gauss curvature calculation
There are many methods for estimating the discrete curvature of triangular meshes. In this paper, the Gauss curvature K of the vertex is calculated on the basis of the triangular mesh model of the skull, based on the method proposed by [11] and so on. The basic idea of this method is to consider a smooth surface as a linear approximation of a family of grids, and to consider the metric properties of each vertex in a triangular mesh as an average measure of the first -order neighbourhood of its voxel model.

Feature Correspondence of Skull
In this paper, the reference skull is represented by T , and the corresponding skull is represented by X . First, the T and X are corrected by Frankfurt coordinate, and the feature points of T are calibrated, then the T and X are registered by TPS, so that the T and X coincide approximately. Finally, the point correspondences between T and X are established based on the constraints of multi geometric features. Specific procedures are as follows: (1) Skull model is voxel oriented. The T and X are voxylated and voxel models are generated by v T and v X respectively. (2) Calculate the normal and curvature. The differential properties of each vertex in v T and v X are calculated: average normal and Gauss curvature, and the descending order of vertices is according to the size of the middle curvature.
(3) Determine the scope of the search. In turn, the vertex of the maximum curvature value in T , which is the most obvious feature, is taken as the current corresponding point, is recorded as v , and the serial number of the voxel unit in three directions of its X , Y , Z is calculated, and the voxel number is used to find the maximum and minimum values of the vertex v in the direction of the voxel, and the voxel order of the vertex v in v X is determined according to the maximum and minimum value. Finally, we determine the search scope of vertex v in v X and use search to represent it.
(4) Determination point correspondence. In accordance with the search range search defined in step (3), expand in v X to find the corresponding point in the current v corresponding to the vertex v T . Searching voxels in search in turn, searching for the corresponding voxel ordinal number in v X , obtaining all the vertices in the voxel ordinal number, and then finding the best corresponding points in these vertices.
From the above 4 steps, we can see that when looking for a physiological correspondence point in the vertex of the T, the search range can be reduced by the body modeling of the skull model. After the convex polyhedron constraint, the candidate points which are not satisfied can be further eliminated. Finally, according to the joint action of the mean normal direction and the Gauss curvature, the best of the corresponding points can be completed match.

Experimental Result Analysis
In this paper, the existing skull sample data set of the Institute of visualization technology of Northwestern University is adopted. In the skull sample database in a random sample of 38 sets of data, the establishment of single face model, as a reference data set corresponding algorithm, the design points corresponding to system reference face and corresponding relation between points corresponding to the dough, the 3D points adding voxel model and the corresponding algorithm of geometric feature constraint.

The results of this method
Experiments are carried out on the proposed algorithm, and the corresponding results of the global point of the skull model and the local points of the skull model are shown in Figure 4 and   Figure 4 and Figure 5 can be seen, this algorithm has achieved good results in the skull model of global and local, the overall results from Figure 4 can be seen without mismatch from Figure 5, partial results can be seen in the same position at both ends of the red line, that the match is correct.

Comparative analysis of different methods
In the skull sample database of 38 sets of data were randomly selected in the same data sample, corresponding method of nearest point [2] based on Euclidean distance and curvature and the method based on the minimum weighted distance corresponding to the method of [12] to establish correspondence between the nose area, using this method, the local points correspond to the results shown in figure 6.  Figure 6 (c) method and the method proposed in this paper are similar, so select the two methods as comparison method this method is that can reflect the accuracy of this method and the rationality. As you can see from Figure 6, the nose area in Figure 6 (a), whether the nose tip or the nose of the nose, meets the physiological consistency of the local point correspondence, the end point of the nasal bone corresponds to the end of the nose, and the alar point corresponds to the alar point. In Figure 6 (b) and Figure 6 (c), some vertices of the nasal alar can be seen clearly and the exact corresponding point is not found. The nasal alar point corresponds to the end point of the nasal bone and mismatches appear.

The accuracy and time consumption of different methods
In order to compare the more intuitive these three methods, 38 sets of data were randomly selected, respectively using the above three methods to establish corresponding relationship, and according to the experimental results, a corresponding number of correct statistics points, calculate the average accuracy rate. The average accuracy is expressed by  Where n represents the number of test samples, N refers to the number of vertices in the reference skull, and i N represents the number of reference points for the reference skull and the first I to match the skull. The accuracy and time consumption of different methods are shown in Table 1.
The average can be seen from the table the method accuracy and time-consuming significantly better than the other two methods, the accuracy of this method is as high as 87.82%, but only time-consuming Euclidean distance method and weighted method to 1/4, curvature method 1/5, so this method has higher accuracy and efficiency. Due to the use of Euclidean distance calculation of all vertices in the corresponding vertex face template to the corresponding wrapper need corresponding nearest point method based on Euclidean distance, and then take the shortest distance as its corresponding points, the calculation process is very time-consuming, greatly increases the time complexity of this algorithm, and there is no screening of the corresponding the vertex, prone to mismatch phenomenon, and the skull model is complex and the corresponding results of this method of unstable; The calculation process of the corresponding to the minimum weighted distance method and method based on curvature is also very time-consuming, and the similarity of local geometric features, geometric constraints and normal curvature is too small, is not conducive to improving the accuracy of the corresponding points, the corresponding results of this method leads to the low accuracy; This paper applies the voxel model, the average method based on Gauss curvature, and convex polyhedron to multi-scale constraint point correspondence algorithm, not only reduce the template for dough corresponding vertex of the corresponding point in the search scope to the corresponding wrapper, improve search efficiency and reduce the time complexity, but also ensure the local area physiological consistent point corresponding to the stability, not because of differences between the model dough and seriously affect the accuracy of the corresponding point. It is proved by experiments that the method is accurate and efficient and has good performance when the skull is the same as the skull to be matched.

Conclusion
A three-dimensional skull correspondence algorithm based on voxel models and multi-geometric feature constraints is proposed in this paper, on the basis of the skull on the deformation of TPS through the skull model voxelization, calculation of vertex differential properties as well as in multi scale constrained convex polyhedron, complete reference is matched with the physiological skull skull consistent point correspondence, and compared with the existing methods. The method of reference point on the skull to find the corresponding point, is in the geometric feature constraints, find the corresponding points in a smaller range, do not need to traverse all the points corresponding to the model, low cost and high efficiency algorithm. The experimental results show that the accuracy of the algorithm is 87.82%, and the time consuming is only 3.15min. The accuracy of point correspondences is improved and the time complexity is reduced.The accuracy of craniofacial restoration is closely related to the accuracy of point correspondences, so the accuracy of point correspondences needs to be further improved in future work; In addition, the weighted distance in normal and curvature, is determined on the basis of experience to determine the weights, and different physiological regions using the same weighting methods, the next step of research so we need to find a method of determining the weights automatically according to different physiological regions.
International Cooperation Project(2013KW04-04) and Shaanxi Provincial Natural Science Basic Research Project(2018JM6061) and the Graduate Scientific Research Foundation of Northwest University (no. YZZ17181).