# A Synthetic Feature Skull Descriptor for 3D Skull Similarity Measurement

Authors: Dan Zhang; Kang Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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 [1–3], craniofacial reconstruction [4–6], 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 [11–14]. 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 [15–18]. 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 [19–21]. 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 λi and eigenfunctions ϕi 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 Figure1. 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 theoreticallyFigure 1
The pipeline of 3D skull similarity measurement.The remainder of this paper is organized as follows. In Section2, 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.


## 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 [11–14]. 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 [15–18]. 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 [19–21]. 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 λi and eigenfunctions ϕi 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 Figure1. 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 theoreticallyFigure 1
The pipeline of 3D skull similarity measurement.The remainder of this paper is organized as follows. In Section2, 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. Figure2 shows some examples from craniofacial database.Figure 2
Some examples of skulls and faces in 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 manifoldM. Let M be a two-dimensional smooth compact manifold with boundary equipped with a Riemann metric, and let M,d 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 C2 real-valued function f in Euclidean space:(1)Δf=∇⋅∇f=∇2f=∂2f∂x2+∂2f∂y2+∂2f∂z2.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, ifg is the metric tensor on M and G is the determinant of the matrix gij, then the LBO can be expressed as [27](2)Δf=∇⋅∇f=1G∑i,j=1ngij∂∂xiGgij∂f∂xj.Since the LBO is self-adjoint and semipositive definite, the LBO onM is decomposed into the matrix product of eigenvalue and eigenfunction: ΔMϕi=λiϕi, where λi is the ith eigenvalue, and ϕi 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 λ2. 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 shapeM 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 vi of the vertex on mesh can be calculated:(3)LBOfvi=∑j=1nωijfvi−fvj.Equation (3) represents a triangular surface sketch of a vertex vi, where fvi and fvj denote the scalar function values defined on M, wij are weights, for vertex vi, and the discrete LBO can be calculated by the method of [27].


## 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. Figure2 shows some examples from craniofacial database.Figure 2
Some examples of skulls and faces in 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 manifoldM. Let M be a two-dimensional smooth compact manifold with boundary equipped with a Riemann metric, and let M,d 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 C2 real-valued function f in Euclidean space:(1)Δf=∇⋅∇f=∇2f=∂2f∂x2+∂2f∂y2+∂2f∂z2.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, ifg is the metric tensor on M and G is the determinant of the matrix gij, then the LBO can be expressed as [27](2)Δf=∇⋅∇f=1G∑i,j=1ngij∂∂xiGgij∂f∂xj.Since the LBO is self-adjoint and semipositive definite, the LBO onM is decomposed into the matrix product of eigenvalue and eigenfunction: ΔMϕi=λiϕi, where λi is the ith eigenvalue, and ϕi 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 λ2. 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 shapeM 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 vi of the vertex on mesh can be calculated:(3)LBOfvi=∑j=1nωijfvi−fvj.Equation (3) represents a triangular surface sketch of a vertex vi, where fvi and fvj denote the scalar function values defined on M, wij are weights, for vertex vi, 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 signatureWKSx,eN of the particle is given by(4)WKSx,eN=Ce∑kϕk2xe−eN−logλk2/2σ2,where x is a point on a manifold, eN is the energy scale parameter, where eN=logE, λk is the kth eigenvalue of LBO, σ is the variance, and Ce=∑ke−eN−logλk2/2σ2−1 is the regularized WKS function.Based on the method of [28], we define the L2 distance of the WKS value between any two points on the shape as the wave kernel distance:(5)dwks2x,y=Ce∑i=1N⟶∞e−eN−logλi2/2σ2ϕy−ϕx2.Assuming that there areNumM 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 nM∗nM symmetric matrix, and the main diagonal element is 0 dwksx,x=0. The wave kernel distance matrix DM is as follows:(6)D=0dwks1,2…dwks1,nMdwks2,10…dwks2,nM…dwksnM,1dwksnM,2…0.If we applyD for measuring the skull similarity, it would cause very high computational complexity. For example, when NumM 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 matrixD and the CDF of D to construct the SWKDD.


### 3.2. The SVD of Wave Kernel Distance Distribution

We first perform singular value decomposition onD to extract the local features of D. Let D∈Rm∗n, and SVD expresses D as a product of the orthonormal matrix U∈Rm∗n and V∈Rm∗n.U is the left singular value matrix of D and V is the right singular value matrix of D; it is well known that(7)D=U∑VT,where UUT=I, VVT=I, and the diagonal elements of ∑ are the nonnegative square roots of the eigenvalues of DTD; they are called singular values. Λ=diag∑=σ1,σ2,…,σn, σ1≥σ2≥⋯≥σr>0,σr+1=⋯=σn. The nonzero singular values in the descending order in this paper are recorded as(8)Λ′=σ1,σ2,…,σr.At the same time, in order to increase the difference in singular values of different skulls, we perform exponential mapping forΛ′.(9)wkdd_svdM=eΛ′.


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

While extracting the features of matrixD 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 Dmin, the maximum value be Dmax, and the interval Dmin,Dmax be divided into I segments. To keep the frequency positive, we shift the coordinate axis to Dmin to eliminate negative values caused by normalization. Therefore, the frequency histogram can be defined as(10)pwksδ=∫χdwksx,y≤δdμxdμy,(11)δ=Dmax−DminI,where pwksδ 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(12)wkdd_cdfM=pwksδμ2M.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:(13)wkddM=wkdd_svdM∗wkdd_cdfM.In order to improve the accuracy and efficiency of calculation, we useZ-score normalization, where μwkddM is the mean value of WKDD of M and σwkddM is the standard deviation of the wkdd(M); wkddM denotes the wkdd(M) after normalization:(14)SWKDDM=wkddM−μwkddMσwkddM.


### 3.4. Skull Similarity Measurement

After obtaining theSWKDDM of the skull, we can measure the similarity between a pair of skulls.(15)dskullp,skullq=∑i=1nSWKDDpi−SWKDDqi,where pi and qi stand for the corresponding point i of SWKDD of skullq and skullq and n is the number of sampling points of SWKDD of skullq and skullq.Thedskullp,skullq distance satisfies the following properties:(i)
Nonnegativity:dskullp,skullq>=0(ii)
Nullity:dskullp,skullq=0 if and only if p=q(iii)
Symmetry:dskullp,skullq = dskullq,skullp(iv)
Triangle inequality:dskullp,skullh+dskullh,skullq>dskullp,skullqIn order to fit the concept of similarity, actually, the more similar the two shapes are, the larger the similarity measure is. We map exponentiallydskullp,skullq and get Dskullp,skullq as the skull similarity measure. The bigger the Dskullp,skullq values are, the more similar the skull p and q are. When q=p, Dskullp,skullq=1, that means skullq and skullq are the same skull.(16)Dskullp,skullq=1EXPdskullp,skullq.The SWKDD for calculation 3D skull similarity algorithm is shown in Algorithm1.Algorithm 1:SWKDD for skull similarity measure algorithm.
Input: skullp, skullq
Ensure: wave kernel distance dwksp,q in equation (6)
Output: distance Dskullp,skullq
(1)
Compute cotan Laplacian (L, B), which defines the Laplacian L′≔B−1L
(2)
ComputeWKSp; WKSq
(3)
Compute the wave kernel distance matrixDp
(4)
Compute wave the kernel distance matrixDq
(5)
Dp⟵Dp−meanDp/stdDp //Normalize the wave kernel distance matrix Dp
(6)
Dq⟵Dq−meanDq/stdDq // Normalize the wave kernel distance matrix DN
(7)
While Dp≠ϕ and Dq≠ϕ do
(8)
Select the largest valueDpmax, Dqmax from Dp and Dq
(9)
Select the smallest valueDpmin, Dqmin from Dp and Dq
(10)
δp=Dpmax−Dpmin/Np; δq=Dqmax−Dqmin/Nq
(11)
Compute the SVD of WKDDwkdd_svdp and wkdd_svdq
(12)
Compute the CDF of WKDDwkdd_cdfp and wkdd_cdfq
(13)
SWKDDp⟵wkdd_svdp∗wkdd_cdfp // computer the SWKDDp
(14)
SWKDDq⟵wkdd_svdq∗wkdd_cdfq // computer the SWKDDq
(15)
Compute the skull similarity distance:Dskullp,skullq
(16)
end while


## 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 signatureWKSx,eN of the particle is given by(4)WKSx,eN=Ce∑kϕk2xe−eN−logλk2/2σ2,where x is a point on a manifold, eN is the energy scale parameter, where eN=logE, λk is the kth eigenvalue of LBO, σ is the variance, and Ce=∑ke−eN−logλk2/2σ2−1 is the regularized WKS function.Based on the method of [28], we define the L2 distance of the WKS value between any two points on the shape as the wave kernel distance:(5)dwks2x,y=Ce∑i=1N⟶∞e−eN−logλi2/2σ2ϕy−ϕx2.Assuming that there areNumM 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 nM∗nM symmetric matrix, and the main diagonal element is 0 dwksx,x=0. The wave kernel distance matrix DM is as follows:(6)D=0dwks1,2…dwks1,nMdwks2,10…dwks2,nM…dwksnM,1dwksnM,2…0.If we applyD for measuring the skull similarity, it would cause very high computational complexity. For example, when NumM 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 matrixD and the CDF of D to construct the SWKDD.


## 3.2. The SVD of Wave Kernel Distance Distribution

We first perform singular value decomposition onD to extract the local features of D. Let D∈Rm∗n, and SVD expresses D as a product of the orthonormal matrix U∈Rm∗n and V∈Rm∗n.U is the left singular value matrix of D and V is the right singular value matrix of D; it is well known that(7)D=U∑VT,where UUT=I, VVT=I, and the diagonal elements of ∑ are the nonnegative square roots of the eigenvalues of DTD; they are called singular values. Λ=diag∑=σ1,σ2,…,σn, σ1≥σ2≥⋯≥σr>0,σr+1=⋯=σn. The nonzero singular values in the descending order in this paper are recorded as(8)Λ′=σ1,σ2,…,σr.At the same time, in order to increase the difference in singular values of different skulls, we perform exponential mapping forΛ′.(9)wkdd_svdM=eΛ′.


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

While extracting the features of matrixD 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 Dmin, the maximum value be Dmax, and the interval Dmin,Dmax be divided into I segments. To keep the frequency positive, we shift the coordinate axis to Dmin to eliminate negative values caused by normalization. Therefore, the frequency histogram can be defined as(10)pwksδ=∫χdwksx,y≤δdμxdμy,(11)δ=Dmax−DminI,where pwksδ 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(12)wkdd_cdfM=pwksδμ2M.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:(13)wkddM=wkdd_svdM∗wkdd_cdfM.In order to improve the accuracy and efficiency of calculation, we useZ-score normalization, where μwkddM is the mean value of WKDD of M and σwkddM is the standard deviation of the wkdd(M); wkddM denotes the wkdd(M) after normalization:(14)SWKDDM=wkddM−μwkddMσwkddM.


## 3.4. Skull Similarity Measurement

After obtaining theSWKDDM of the skull, we can measure the similarity between a pair of skulls.(15)dskullp,skullq=∑i=1nSWKDDpi−SWKDDqi,where pi and qi stand for the corresponding point i of SWKDD of skullq and skullq and n is the number of sampling points of SWKDD of skullq and skullq.Thedskullp,skullq distance satisfies the following properties:(i)
Nonnegativity:dskullp,skullq>=0(ii)
Nullity:dskullp,skullq=0 if and only if p=q(iii)
Symmetry:dskullp,skullq = dskullq,skullp(iv)
Triangle inequality:dskullp,skullh+dskullh,skullq>dskullp,skullqIn order to fit the concept of similarity, actually, the more similar the two shapes are, the larger the similarity measure is. We map exponentiallydskullp,skullq and get Dskullp,skullq as the skull similarity measure. The bigger the Dskullp,skullq values are, the more similar the skull p and q are. When q=p, Dskullp,skullq=1, that means skullq and skullq are the same skull.(16)Dskullp,skullq=1EXPdskullp,skullq.The SWKDD for calculation 3D skull similarity algorithm is shown in Algorithm1.Algorithm 1:SWKDD for skull similarity measure algorithm.
Input: skullp, skullq
Ensure: wave kernel distance dwksp,q in equation (6)
Output: distance Dskullp,skullq
(1)
Compute cotan Laplacian (L, B), which defines the Laplacian L′≔B−1L
(2)
ComputeWKSp; WKSq
(3)
Compute the wave kernel distance matrixDp
(4)
Compute wave the kernel distance matrixDq
(5)
Dp⟵Dp−meanDp/stdDp //Normalize the wave kernel distance matrix Dp
(6)
Dq⟵Dq−meanDq/stdDq // Normalize the wave kernel distance matrix DN
(7)
While Dp≠ϕ and Dq≠ϕ do
(8)
Select the largest valueDpmax, Dqmax from Dp and Dq
(9)
Select the smallest valueDpmin, Dqmin from Dp and Dq
(10)
δp=Dpmax−Dpmin/Np; δq=Dqmax−Dqmin/Nq
(11)
Compute the SVD of WKDDwkdd_svdp and wkdd_svdq
(12)
Compute the CDF of WKDDwkdd_cdfp and wkdd_cdfq
(13)
SWKDDp⟵wkdd_svdp∗wkdd_cdfp // computer the SWKDDp
(14)
SWKDDq⟵wkdd_svdq∗wkdd_cdfq // computer the SWKDDq
(15)
Compute the skull similarity distance:Dskullp,skullq
(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 isOn3, 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 wkdd_cdf, wkdd_svd, 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 theSWKDD. We randomly select 10 skulls of males from the database. Figure 3 shows the wkdd_cdf and SWKDD of male skulls. From Figure 3, we can see that SWKDD can better distinguish between different skulls compared with wkdd_cdf. Because SWKDD 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.Figure 3
Thewkdd_cdf and SWKDD of male skulls. (a) Male skulls’ wkdd_cdf value. (b) Male skulls’ SWKDD value.
(a)
(b)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 Section4.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 R1 and R2, respectively. Based on the line interpolation deformation process between the two models, the middle four models of the R2 model generated from the R1 deformation are selected.We deform the modelR1 and R2 (for both male and female) simulating and then produce morph-deformed skulls and faces using the formula Fi=1−a∗R1+a∗R2 a=0.2,0.4,0.6,0.8, where the parameter a can control the degree upon how much the two skulls or faces are mixed: if a=0, the generated face is the same with F1. If a=1, it is the same with F2. The value of a reflects the similarity between the generated model with the original models. Take a=0.2, for example, the generated face should be closer to the skull F1, while the face generated by a=0.8 should be closer to the face F2. F1=0.8∗R1+0.2∗R2, F1=0.6∗R1+0.4∗R2, F1=0.4∗R1+0.6∗R2, and M1=0.2∗R1+0.8∗R2.R
1 and R2 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.Figure 4
Similarity results of males and females.Tables1 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.Table 1
The skull similarity between every two morph female skulls byR1 and R2.
R1
F1
F2
F3
F4
R2
R1
1.000000
0.819253
0.711221
0.697863
0.695526
0.677464
F1
0.819253
1.000000
0.973012
0.951868
0.953956
0.932365
F2
0.711221
0.973012
1.000000
0.991893
0.991472
0.980058
F3
0.697863
0.951868
0.991893
1.000000
0.999588
0.997192
F4
0.695526
0.953956
0.991472
0.999588
1.000000
0.996961
R2
0.677464
0.932365
0.980058
0.997192
0.996961
1.000000Table 2
The skull similarity between every two morph male skulls byR1 and R2.
R1
F1
F2
F3
F4
R2
R1
1.000000
0.998806
0.997007
0.991150
0.973169
0.819047
F1
0.998806
1.000000
0.998545
0.988778
0.967815
0.826089
F2
0.997007
0.998545
1.000000
0.982464
0.957577
0.842712
F3
0.991150
0.988778
0.982464
1.000000
0.993286
0.748852
F4
0.973169
0.967815
0.957577
0.993286
1.000000
0.690204
R2
0.819047
0.826089
0.842712
0.748852
0.690204
1.000000Figure 5
Similarity results of males and females.


### 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.Figure6 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].Figure 6
Similarity results of males and females. (a) The skull similarity and face similarity of male. (b) The skull similarity and face similarity of female.
(a)
(b)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 Figures7 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.Figure 7
Skull similarity and face similarity thermogram results of males.Figure 8
Skull similarity and face similarity thermogram results of females.


### 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. Figure9 shows the different sampling of a female skull.Figure 9
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 Table3. 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.Table 3
The skull similarity between every two resample female skulls.
Skull1
Resample 50%
Resample 25%
Resample 12.5%
Resample 6.5%
Skull1
1.000000
0.9587798
0.9056897
0.861239
0.813233
Resample 50%
0.9587798
1.000000
0.938345
0.888456
0.847811
Resample 25%
0.9056897
0.938345
1.000000
0.922565
0.908961
Resample 12.5%
0.861239
0.888456
0.922565
1.000000
0.981453
Resample 6.5%
0.813233
0.847811
0.908961
0.981453
1.000000


### 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 usingwkdd_cdf, wkdd_svd, 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. Figure10 shows the ROC curves of the gender classification both of the 3D face and skull using wkdd_cdf, wkdd_svd, and SWKDD separately. From Figure 10, the AUC values of SWKDD are bigger than those of wkdd_cdf and wkdd_svd both of the 3D face and skull, and these mean that SWKDD has good performance than wkdd_cdf and wkdd_svd under unsupervised gender classification.Figure 10
The ROC curves of the gender classification both of the 3D face and skull.


## 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 R1 and R2, respectively. Based on the line interpolation deformation process between the two models, the middle four models of the R2 model generated from the R1 deformation are selected.We deform the modelR1 and R2 (for both male and female) simulating and then produce morph-deformed skulls and faces using the formula Fi=1−a∗R1+a∗R2 a=0.2,0.4,0.6,0.8, where the parameter a can control the degree upon how much the two skulls or faces are mixed: if a=0, the generated face is the same with F1. If a=1, it is the same with F2. The value of a reflects the similarity between the generated model with the original models. Take a=0.2, for example, the generated face should be closer to the skull F1, while the face generated by a=0.8 should be closer to the face F2. F1=0.8∗R1+0.2∗R2, F1=0.6∗R1+0.4∗R2, F1=0.4∗R1+0.6∗R2, and M1=0.2∗R1+0.8∗R2.R
1 and R2 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.Figure 4
Similarity results of males and females.Tables1 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.Table 1
The skull similarity between every two morph female skulls byR1 and R2.
R1
F1
F2
F3
F4
R2
R1
1.000000
0.819253
0.711221
0.697863
0.695526
0.677464
F1
0.819253
1.000000
0.973012
0.951868
0.953956
0.932365
F2
0.711221
0.973012
1.000000
0.991893
0.991472
0.980058
F3
0.697863
0.951868
0.991893
1.000000
0.999588
0.997192
F4
0.695526
0.953956
0.991472
0.999588
1.000000
0.996961
R2
0.677464
0.932365
0.980058
0.997192
0.996961
1.000000Table 2
The skull similarity between every two morph male skulls byR1 and R2.
R1
F1
F2
F3
F4
R2
R1
1.000000
0.998806
0.997007
0.991150
0.973169
0.819047
F1
0.998806
1.000000
0.998545
0.988778
0.967815
0.826089
F2
0.997007
0.998545
1.000000
0.982464
0.957577
0.842712
F3
0.991150
0.988778
0.982464
1.000000
0.993286
0.748852
F4
0.973169
0.967815
0.957577
0.993286
1.000000
0.690204
R2
0.819047
0.826089
0.842712
0.748852
0.690204
1.000000Figure 5
Similarity results of males and females.


## 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.Figure6 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].Figure 6
Similarity results of males and females. (a) The skull similarity and face similarity of male. (b) The skull similarity and face similarity of female.
(a)
(b)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 Figures7 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.Figure 7
Skull similarity and face similarity thermogram results of males.Figure 8
Skull similarity and face similarity thermogram results of females.


## 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. Figure9 shows the different sampling of a female skull.Figure 9
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 Table3. 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.Table 3
The skull similarity between every two resample female skulls.
Skull1
Resample 50%
Resample 25%
Resample 12.5%
Resample 6.5%
Skull1
1.000000
0.9587798
0.9056897
0.861239
0.813233
Resample 50%
0.9587798
1.000000
0.938345
0.888456
0.847811
Resample 25%
0.9056897
0.938345
1.000000
0.922565
0.908961
Resample 12.5%
0.861239
0.888456
0.922565
1.000000
0.981453
Resample 6.5%
0.813233
0.847811
0.908961
0.981453
1.000000


## 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 usingwkdd_cdf, wkdd_svd, 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. Figure10 shows the ROC curves of the gender classification both of the 3D face and skull using wkdd_cdf, wkdd_svd, and SWKDD separately. From Figure 10, the AUC values of SWKDD are bigger than those of wkdd_cdf and wkdd_svd both of the 3D face and skull, and these mean that SWKDD has good performance than wkdd_cdf and wkdd_svd under unsupervised gender classification.Figure 10
The ROC curves of the gender classification both of the 3D face and skull.


## 5. Discussion

The abovementioned experimental results show the efficiency of our proposed skull similarity measurement based on SWKDD. From Figures3 and 10, compared with wkdd_cdf and wkdd_svd, the SWKDD 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 SWKDD 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 5–8, 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.



---

*Source: 8083504-2019-11-30.xml*