Abstract

3D skull similarity measurement is a challenging and meaningful task in the fields of archaeology, forensic science, and anthropology. However, it is difficult to correctly and directly measure the similarity between 3D skulls which are geometric models with multiple border holes and complex topologies. In this paper, based on the synthetic feature method, we propose a novel 3D skull descriptor, synthetic wave kernel distance distribution (SWKDD) constructed by the laplace–beltrami operator. By defining SWKDD, we obtain a concise global skull representation method and transform the complex 3D skull similarity measurement into a simple 1D vector similarity measurement. First, we give the definition and calculation of SWKDD and analyse its properties. Second, we represent a framework for 3D skull similarity measurement using the SWKDD of 3D skulls and details of the calculation steps involved. Finally, we validate the effectiveness of our proposed method by calculating the similarity measurement of 3D skulls based on the real craniofacial database.

1. Introduction

Recently, with the continuous development of craniofacial morphology and biological computing, the research of 3D skull has attracted the attention of researchers. The 3D skull is an intrinsic creature of the human face, which is a well preserved bone under the effect of fire, humidity, and temperature changes. It can be widely applied in the fields of craniofacial identification [13], craniofacial reconstruction [46], and so on. The similarity measurement of the skull is one of the key techniques in the abovementioned craniofacial research. It can be used to determine the similarity interval of craniofacial identification or evaluation of craniofacial reconstruction.

1.1. Related Works

At present, researchers have not done enough research on 3D skull similarity measurement, especially directly measuring methods. The main reasons is as follows: first, it is not easy for researchers to acquire and process 3D skull models; second, the topological structure of the 3D skull is more complex than that of the 3D face; and finally, because the shape of different 3D skull models are very similar, although human visual system can quickly judge the similarity of two 3D faces in the brain, it is difficult for the visual system to automatically judge the similarity of two 3D skulls.

Delfino et al. [7] used an original shape analytical morphometry (SAM) software package to match a skull and a photograph of a face. Zhu and Geng [8] proposed an approach of 3D craniofacial geometry shape similarity estimation based on principal warps. Zhao et al. [9] used the sparse principal component analysis (SPCA) method to evaluate the similarity between two sets of craniofacial models. Zhang et al. [10] proposed a local shape descriptor Harmonic Wave Kernel Signature (HWKS) to measure 3D skull similarity by calculating the HWKS of the corresponding points of a pair of 3D skulls. Some works transformed the similarity of skulls into the similarity of its corresponding reconstructed 3D faces through craniofacial reconstruction [1114]. The more similar the corresponding reconstructed faces are, the more similar the skulls are, but this similarity result is very dependent on the craniofacial reconstruction method. These abovementioned methods cannot directly measure the similarity between the skulls, and the results obtained are not accurate.

Some works are based on the statistical model methods to extract the features of skulls and faces [1518]. Duan et al. [17] proposed a skull identification method based on PCA that matches an unknown skull with enrolled 3D faces, in which the mapping between the skull and face is obtained using canonical correlation analysis. Shui et al. [18] investigated the extent to which the choice of principal components (PCs) affects the analysis of the craniofacial relationships and sexual dimorphism. In contrast, in practical applications, the principal components with small contributions also contain explicit information that distinguishes different skulls.

Skull as a special complex 3D shape, researchers used 3D shape descriptors to extract the characteristics of skulls [1921]. Jin et al. [19] used canny operator and the sliding window method for edge detection and boundary tracking to measure the skull similarity measurement. Yu et al. [20] used the heat kernel signature to express and repair the 3D digital skull model. The heat kernel signature (HKS) [21] is derived from the eigenvalues and eigenfunctions of the Laplace–Beltrami operator (LBO) [22] on the surface of a Riemann manifold. Liu [23] proposed a method to measure 3D skull similarity by calculating the Fréchet distance of the radial curve on the skull with filling the multiple border holes of skulls. However, this method is not very accurate because the hole filling operation ignores the details of border holes of skulls.

1.2. Contribution

In this paper, we propose a global shape descriptor, synthetic wave kernel distance distribution (SWKDD), to construct the 3D skull similarity measurement. First, our method can directly measure the similarity between a pair of skulls without filling multiple border holes. We choose the wave kernel signature (WKS) to define a concise 3D skull descriptor. The WKS [24] uses a bandpass filter to clearly separate different sets of frequencies on the shape and well describe the complex skull model. Second, our method can avoid finding the same number of the corresponding points on skulls before measuring the skull similarity. We calculate the cumulative distribution function and singular value of distance of the WKS based on the idea of Brostein and Bronstein [25]. Bronstein applied statistical methods to select shape diffusion distance and commute-time distance, calculated the distribution histogram of the diffusion distance and commute-time distance, and then used the histogram as a global shape descriptor.

We construct a pipeline of 3D skull similarity measurement which schematically describes the generic framework for measuring the similarity of 3D skulls based on the SWKDD, as shown in Figure 1. In fact, this framework also can be applied to 3D face similarity measurement. In summary, the contributions of our study are as follows:(i)A novel shape descriptor SWKDD involved in LBO is defined, which is an intrinsic descriptor, independent of the coordinate and without filling multiple border holes of skulls(ii)A framework for directly measuring the similarity of the 3D skull is proposed, and we map the 3D skull into the 1D vector space which is constructed by the SWKDD, which is convenient to measure the similarity between a pair of skulls(iii)The effectiveness of our method is validated by different experimental strategies on the real Asian Mongolian craniofacial database, and we also analyze the experimental results theoretically

The remainder of this paper is organized as follows. In Section 2, we introduce the materials and fundamentals for our method. In Section 3, we provide the SWKDD in detail and construct a 3D skull similarity measurement. In Section 4, we show our experimental results. Finally, we draw conclusions regarding our study in Section 5.

2. Materials and Fundamentals

2.1. Materials

The Asian Mongolian craniofacial database in this article was obtained from 140 volunteers (70 male and 70 female) who mostly came from the Han ethnic group of Northern China in Xianyang Hospital located in western China [6, 26]. The database was established by scanning the heads with a clinical multislice CT scan system (Siemens Sensation16). When the data were collected, the volunteers were supine, and their hands were naturally sagging. Their feet were close together, and their eyes were facing forward. The ages of volunteers ranged from 19 to 75. Axial position helical scanning was used to reconstruct a thickness of 0.75 mm. The collected data were stored in Standard DICOM 3.0 images with an interlayer slice resolution of 512 ∗ 512 for all samples. The color depth is 16 bits. The collected data is stored according to DICOM 3.0 standard. One scan requires about 320 CT images to obtain complete head image data.

Similarity measurements of 3D faces and skulls are usually aimed at a complete 3D human face or skull model. After filling tiny holes, removing scattered points and nonrigid registration, the 3D craniofacial model has a complete manifold structure. These models have a unified coordinate system and data scale to eliminate the impact of data size and posture. Figure 2 shows some examples from craniofacial database.

In this section, we introduce the craniofacial materials and fundamentals for 3D skull similarity measurement, and our method is based on the LBO. The LBO is a well-known intrinsic operator which can be decomposed by spectral decomposition. In the spectral analysis method, the eigenfunction and eigenvalue of the LBO can be used in defining different spectral shape descriptors and distances. By calculating spectral distance distributions of a pair of shapes, the local or global matching of a pair of shapes can be compared, and a similarity result between a pair of 3D skulls is obtained.

2.2. LBO

To effectively represent the intrinsic information and geometric features of the shape, we consider the 3D shape as a manifold M. Let M be a two-dimensional smooth compact manifold with boundary equipped with a Riemann metric, and let be a metric space. For a compact Riemann manifold M, we apply the spectral distance method to define shape features to construct feature space, which is in turn closely connected with the notion of LBO. The Laplace operator is a differential operator defined by the gradient and divergence of a real-valued function f in Euclidean space:

By equipping the Laplace operator with a Riemann manifold metric, we obtain the LBO. According to the definition of the Riemann manifold gradient and divergence, if is the metric tensor on M and G is the determinant of the matrix , then the LBO can be expressed as [27]

Since the LBO is self-adjoint and semipositive definite, the LBO on M is decomposed into the matrix product of eigenvalue and eigenfunction: , where is the eigenvalue, and is the corresponding eigenfunction. If the Neumann boundary condition is used in a closed region, then the first eigenvalue is 0, and the smallest nonzero eigenvalue is . The LBO can be analytically calculated for some geometrical shape (e.g., rectangular or cylindrical). The LBO eigenfunctions are intrinsic to the manifold, and the ones related to smaller eigenvalues correspond to smooth and slowly varying functions.

For numerical computation, the shape M can be represented by a finite set of points. In discrete mathematics, the finite-dimensional discrete LBO is typically called the discrete Laplace–Beltrami matrix. On a triangular mesh with a vertex number of n, the discrete LBO at the of the vertex on mesh can be calculated:

Equation (3) represents a triangular surface sketch of a vertex , where and denote the scalar function values defined on M, are weights, for vertex , and the discrete LBO can be calculated by the method of [27].

3. Wave Kernel Distance Distribution and 3D Skull Similarity Measurement

In this section, we construct a global shape descriptor, synthetic wave kernel distance distribution (SWKDD), and propose a 3D skull similarity measurement based on SWKDD. We calculate the cumulative distribution function and singular value of the distance of the WKS to construct the SWKDD.

3.1. Wave Kernel Distance Distribution

For each point on a shape, a shape descriptor called the wave kernel signature (WKS) operator is defined by measuring the average probability distribution of quantum particles with different energy levels. The wave kernel signature of the particle is given bywhere x is a point on a manifold, is the energy scale parameter, where , is the eigenvalue of LBO, σ is the variance, and is the regularized WKS function.

Based on the method of [28], we define the distance of the WKS value between any two points on the shape as the wave kernel distance:

Assuming that there are vertices on the manifold M (vertex number starting from 1), if we seek the wave kernel distance for any two points on M, the wave kernel distance matrix of the shape is an symmetric matrix, and the main diagonal element is 0 . The wave kernel distance matrix is as follows:

If we apply D for measuring the skull similarity, it would cause very high computational complexity. For example, when is 10000, D is a very large high-dimensional matrix. So, when we measure the skull similarity, we need to reduce the dimensionality of D.

We consider constructing a novel shape descriptor based on statistical theory and matrix decomposition. The singular value decomposition (SVD) maps variables into low-dimensional spaces and preserves important and efficient local features of variables. The main idea of SVD is to decompose the matrix into several matrix multiplications, which is an effective method for extracting features. At the same time, the cumulative distribution function (CDF) can fully describe the probability distribution of a random variable and the global statistical properties of variables.

We hope that the proposed descriptor can preserve the geometric and statistical properties of the skull as much as possible. Therefore, we calculate the SVD of the wave kernel distance matrix D and the CDF of D to construct the SWKDD.

3.2. The SVD of Wave Kernel Distance Distribution

We first perform singular value decomposition on D to extract the local features of D. Let , and SVD expresses D as a product of the orthonormal matrix and .

U is the left singular value matrix of D and V is the right singular value matrix of D; it is well known thatwhere , , and the diagonal elements of are the nonnegative square roots of the eigenvalues of ; they are called singular values. , . The nonzero singular values in the descending order in this paper are recorded as

At the same time, in order to increase the difference in singular values of different skulls, we perform exponential mapping for .

3.3. The CDF of Wave Kernel Distance Distribution and the SWKDD

While extracting the features of matrix D by SVD, the statistical features of matrix D is represented. δ is the distance threshold (the threshold is divided by the difference of the distance and the number of frequency histograms, as given in equation (11)), μ is a norm metric defined in M, and χ is the indicator function. Let the minimum value of D be , the maximum value be , and the interval be divided into I segments. To keep the frequency positive, we shift the coordinate axis to to eliminate negative values caused by normalization. Therefore, the frequency histogram can be defined aswhere defined this way is the measure of pairs of points with their distance no larger than δ. In probability theory and statistics, the cumulative distribution function (CDF) is an integral of the probability density function (PDF), which is used to fully specify the probability distribution of random variables. For a manifold M, the CDF of D is calculated by

In machine learning, synthetic feature can improve the discrimination ability of a single feature. A more effective and common strategy is to multiply the values of the two features to form a feature combination [29]. In this paper, we synthesize the SVD and the CDF of D of the skull M as a novel skull descriptor, synthetic wave kernel distance distribution:

In order to improve the accuracy and efficiency of calculation, we use Z-score normalization, where is the mean value of WKDD of M and is the standard deviation of the wkdd(M); denotes the wkdd(M) after normalization:

3.4. Skull Similarity Measurement

After obtaining the of the skull, we can measure the similarity between a pair of skulls.where and stand for the corresponding point i of SWKDD of and and n is the number of sampling points of SWKDD of and .

The distance satisfies the following properties:(i)Nonnegativity: (ii)Nullity: if and only if (iii)Symmetry:  = (iv)Triangle inequality:

In order to fit the concept of similarity, actually, the more similar the two shapes are, the larger the similarity measure is. We map exponentially and get as the skull similarity measure. The bigger the values are, the more similar the skull p and q are. When , , that means and are the same skull.

The SWKDD for calculation 3D skull similarity algorithm is shown in Algorithm 1.

Input: ,
Ensure: wave kernel distance in equation (6)
Output: distance
(1) Compute cotan Laplacian (L, B), which defines the Laplacian
(2) Compute ;
(3) Compute the wave kernel distance matrix
(4) Compute wave the kernel distance matrix
(5) //Normalize the wave kernel distance matrix
(6) // Normalize the wave kernel distance matrix
(7)While and do
(8)  Select the largest value , from and
(9)  Select the smallest value , from and
(10)  ;
(11)  Compute the SVD of WKDD and
(12)  Compute the CDF of WKDD and
(13)   // computer the
(14)   // computer the
(15)  Compute the skull similarity distance:
(16)end while

4. Experiment Results

We perform experiments on MATLAB 2015 of a 64 bit 32G memory, Win10 system. The time complexity of calculating craniofacial similarity with SWKDD is , and n is the number of sampling points on the 3D craniofacial model. The experiments are constructed by five parts. Firstly, we evaluate the skull similarity measurement experimental verification on morph Data. Secondly, we evaluate the skull similarity measurement experimental verification on real data. Thirdly, we evaluate the resample robustness of our skull similarity measurement. Fourthly, we evaluate the gender classification both of 3D face and skull based on , , and SWKDD separately using support vector machine(SVM) [30]. Finally, we give the discussion of the abovementioned experiments.

Before verifying the skull similarity proposed in this paper, we give an illustration of the effectiveness of the . We randomly select 10 skulls of males from the database. Figure 3 shows the and of male skulls. From Figure 3, we can see that can better distinguish between different skulls compared with . Because introduces the singular value of the cumulative distribution distance of the wave kernel signature, it enhances the characterization of the same type of 3D shapes.

For 3D skulls, people cannot subjectively identify skull similarity. Only the experts who study archaeology and forensic science can distinguish the difference between the skulls. So, we have the following two experiments to verify our method.

First, we have experiments on morph data in Section 4.1; we generate several morph-deformed male and female skulls which we know the similarity between each two skulls in advance. Therefore, we can verify whether the similarity calculated by our method is correct or not. Second, we have experiments on different real craniofacial models and apply our method to calculate whether the trend of skull similarity is consistent with the trend of face similarity both of male and female in Section 4.2.

4.1. The Skull Similarity Measurement Experimental Verification on Morph Data

As in [31], two 3D craniofacial models are randomly selected from the craniofacial database, namely and , respectively. Based on the line interpolation deformation process between the two models, the middle four models of the model generated from the deformation are selected.

We deform the model and (for both male and female) simulating and then produce morph-deformed skulls and faces using the formula , where the parameter a can control the degree upon how much the two skulls or faces are mixed: if , the generated face is the same with . If , it is the same with . The value of a reflects the similarity between the generated model with the original models. Take , for example, the generated face should be closer to the skull , while the face generated by should be closer to the face . , , , and .

and deformation processes are shown in the Figure 4. Since it is difficult to identify the degree of similarity between the skulls subjectively, we give the corresponding face model to facilitate the indirect judgement of its similarity.

Tables 1 and 2 show the similarity of the male skull and female skull, and the matrix is symmetric. The diagonal value of the matrix is 1, indicating the greatest similarity between the same skulls. As can be seen from Tables 1 and 2 , the trend of similarity value between the female skull R1, F1, F2, F3, F4, and R2 with the skull R1 defined is sequentially decreased (also can be seen in Figure 5). The similarity of male skull has the same laws, and the results by morph deformation skull are also consistent with the person’s subjective judgment. In this paper, our proposed SWKDD can represent the local and global properties of the skull; therefore, the 3D skull similarity can be measured correctly and efficiently.

4.2. The Skull Similarity Measurement Experimental Verification on Real Data

In order to verify the skull similarity proposed in our paper, we, respectively, give a line diagram of the skull and face similarity of male and female.

Figure 6 shows a line graph of the similarity change trend both of the males and female’s skulls and faces based on SWKDD (1st to others). When the similarity between the No. 1 skull and other skulls increases, the similarity between the No. 1 face and other faces also increases. From Figure 6, the similarity between the skulls and the similarity between the faces show basically consistent trends for both male and female, indicating the effectiveness of our proposed method. Importantly, we discover the divergence of a male’s skull and face similarity trend is smaller than that of a female. Because a male’s face features are more prominent, and their skull shape can better determine the face shape [32].

In order to display the laws of data more intuitively, we randomly selected 10 groups of craniofacial data for male and female (including 5 groups of male and female), calculated the similarity between the skull and the similarity between the faces both of male and female, and compared the differences between the similarity of the skull and face using the thermodynamic in Figures 7 and 8. Where blue indicates the maximum similarity value and yellow indicates the minimum similarity value, through analysis of the results of the overall change trends in skull similarity and face similarity of males and females, the skull similarity and face similarity are generally related, and the colour change trends are similar, indicating the effectiveness of our proposed method.

4.3. The Resample Robustness of SWKDD

In engineering applications, although the 3D models with high sampling vertex look smoother and the calculated results will be more accurate, it is undeniable that its computational complexity is higher than that of the rough model, which consumes more time and space. Therefore, in practical applications, the algorithm can be applied to models with different numbers of sampling vertex. The 3D skull descriptor SWKDD is mainly constructed by the LBO, which is a smooth second-order partial differential operator and intrinsic properties of manifolds. Therefore, the different sampling of craniofacial models has little effect on our proposed skull similarity measurement, so we can use it to calculate the similarity of 3D models with different numbers of sampling vertex. Figure 9 shows the different sampling of a female skull.

We calculate the similarity between the 3D skull model with different numbers of sampling vertex and the original model by our method, and the results are shown in Table 3. As can be seen from Table 3, as the number of the skull’s sampling vertex decreases, the similarity between the 3D skull model and the original model decreases, but our method can ensure that the similarity between the simplified model and the original model is more than 0.8.

4.4. Unsupervised Gender Classification Based on SWKDD

In order to further show the effectiveness of SWKDD, we use SVM to evaluate the unsupervised gender classification both of 3D face and skull only by using , , and SWKDD separately without labels. An SVM is a discriminative classifier formally defined by a separating hyperplane. It can be used for classification, regression, or other tasks. In this paper, we use LIBSVM tools to conduct our experiments and show the receiver operating characteristic (ROC) curves of the results.

Gender classification is a dichotomy problem, in which cases are divided into positive or negative categories. Figure 10 shows the ROC curves of the gender classification both of the 3D face and skull using , , and SWKDD separately. From Figure 10, the AUC values of SWKDD are bigger than those of and both of the 3D face and skull, and these mean that SWKDD has good performance than and under unsupervised gender classification.

5. Discussion

The abovementioned experimental results show the efficiency of our proposed skull similarity measurement based on SWKDD. From Figures 3 and 10, compared with and , the synthetically constructed by the CDF and SVD of wave kernel signature distance can better describe the global and local feature of the 3D skull and face. Therefore, the can capture and show more differences in 3D skulls with similar contours and topological structure. Based on the SWKDD, we define a 3D skull similarity and can also be used for 3D faces. From Figures 58, we can clearly see that both for male and female, AAata. Moreover, compared with female, the change trend of the skull similarity and face similarity of male is more consistent, indicating that the relationship between the skull and face of male is closer.

6. Conclusion

In this article, we define an efficient shape descriptor SWKDD. By defining the SWKDD, we map the 3D skulls into the 1D vector which characterizes the statistical and geometric properties of 3D skulls. The skull feature extraction and representation are convenient by SWKDD. Also, we provide a 3D skull similarity measurement framework based on the SWKDD. Our framework is more ubiquitous for measuring the similarity of 3D skulls without filling the multiple border holes and finding the landmark point on the 3D skull, and our similarity measurement framework can be used to measure 3D models with different numbers of sampling vertex efficiently. It can be seen from experiments that our proposed measurement can effectively and accurately measure the similarity of the 3D skull and face with different numbers of sampling vertex.

In future work, we will find a solution to solve the problem such as using SWKDD for automatic skull identification and evaluating the craniofacial reconstruction. We will attempt to find more convenient skull features to construct skull representation which includes the global and local characteristics.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by Beijing Natural Science Foundation (194072), Beijing Natural Science Foundation (4194072), Beijing Excellent Talents Training Project (2017000020124G057), and Beijing Municipal Education Commission Science and Technology General Project (KM201910028001).