Detecting unique phenotypes across populations can be achieved by quantitative analysis of bone shape, provided that the databases of normal and abnormal pathologies are available. The recent surge of interest in the spectral analysis of the Laplace-Beltrami operator (LBO) has resulted in a considerable number of spectral shape signatures that have been successfully applied to a broad range of areas, including manifold learning [Belkin:06], object recognition and deformable shape analysis [Levy:06, Reuter:06, Rustamov:07, Bronstein:11, Masoumi:16, Rodola:SHREC17], medical imaging [Chaudhari:14], multimedia protection [Tarmissi:09], and shape classification [Masoumi:17]
. The diversified nature of these applications is a powerful testimony of the practical usage of spectral shapes signatures, which are usually defined as feature vectors representing local and/or global characteristics of a shape and may be broadly classified into two main categories: local and global descriptors. Local descriptors (also called point signatures) are defined on each point of the shape and often represent the local structure of the shape around that point, while global descriptors are usually defined on the entire shape.
Most point signatures may easily be aggregated to form global descriptors by integrating over the entire shape. Rustamov [Rustamov:07]
proposed a local feature descriptor referred to as the global point signature (GPS), which is a vector whose components are scaled eigenfunctions of the LBO evaluated at each surface point. The GPS signature is invariant under isometric deformations of the shape, but it suffers from the problem of eigenfunctions’ switching whenever the associated eigenvalues are close to each other. This problem was lately well handled by the heat kernel signature (HKS)[Sun:09], which is a temporal descriptor defined as an exponentially-weighted combination of the LBO eigenfunctions. HKS is a local shape descriptor that has a number of desirable properties, including robustness to small perturbations of the shape, efficiency and invariance to isometric transformations. The idea of HKS was also independently proposed by Gȩbal et al. [Gebal:09] for 3D shape skeletonization and segmentation under the name of auto diffusion function. From the graph Fourier perspective, it can be seen that HKS is highly dominated by information from low frequencies, which correspond to macroscopic properties of a shape. To give rise to substantially more accurate matching than HKS, the wave kernel signature (WKS) [Aubry:11]
was proposed as an alternative in an effort to allow access to high-frequency information. Using the Fourier transform’s magnitude, Kokkinoset al. [Kokkinos:10] introduced the scale invariant heat kernel signature (SIHKS), which is constructed based on a logarithmically sampled scale-space.
One of the simplest spectral shape signatures is Shape-DNA [Reuter:06], which is an isometry-invariant global descriptor defined as a truncated sequence of the LBO eigenvalues arranged in increasing order of magnitude. Gao et al. [Gao:14] developed a variant of Shape-DNA, referred to as compact Shape-DNA (cShape-DNA), which is an isometry-invariant signature resulting from applying the discrete Fourier transform to the area-normalized eigenvalues of the LBO. Chaudhari et al. [Chaudhari:14] presented a slightly modified version of the GPS signature by setting the LBO eigenfunctions to unity. This signature, called GPS embedding, is defined as a truncated sequence of inverse square roots of the area-normalized eigenvalues of the LBO. A comprehensive list of spectral descriptors can be found in [Lian:13].
On the other hand, wavelet analysis has some major advantages over Fourier transform, which makes it an interesting alternative for many applications. In particular, unlike the Fourier transform, wavelet analysis is able to perform local analysis and also makes it possible to perform a multiresolution analysis. Classical wavelets are constructed by translating and scaling a mother wavelet, which is used to generate a set of functions through the scaling and translation operations. The wavelet transform coefficients are then obtained by taking the inner product of the input function with the translated and scaled waveforms. The application of wavelets to graphs (or triangle meshes in geometry processing) is, however, problematic and not straightforward due in part to the fact that it is unclear how to apply the scaling operation on a signal (or function) defined on the mesh vertices. To tackle this problem, Coifman et al. [Coifman:06] introduced the diffusion wavelets, which generalize the classical wavelets by allowing for multiscale analysis on graphs. The construction of diffusion wavelets interacts with the underlying graph through repeated applications of a diffusion operator, which induces a scaling process. Hammond et al. [Hammond:11] showed that the wavelet transform can be performed in the graph Fourier domain, and proposed a spectral graph wavelet transform that is defined in terms of the eigensystem of the graph Laplacian matrix. More recently, a spectral graph wavelet signature (SGWS) was introduced in [Chunyuan:13b], and it has shown superior performance over HKS and WKS in 3D shape retrieval. SGWS is a multiresolution local descriptor that is not only isometric invariant, but also compact, easy to compute and combines the advantages of both band-pass and low-pass filters.
A popular approach for transforming local descriptors into global representations that can be used for surface analysis is the bag-of-features (BoF) paradigm [Bronstein:11]
. The BoF model represents each object in the dataset as a collection of unordered feature descriptors extracted from local areas of the shape, just as words are local features of a document. A baseline BoF approach quantizes each local descriptor to its nearest cluster center using K-means clustering and then encodes each shape as a histogram over cluster centers by counting the number of assignments per cluster. These cluster centers form a visual vocabulary or codebook whose elements are often referred to as visual words or codewords. Although the BoF paradigm has been shown to provide significant levels of performance, it does not, however, take into consideration the spatial relations between features, which may have an adverse effect not only on its descriptive ability but also on its discriminative power. To account for the spatial relations between features, Bronsteinet al. introduced a generalization of a bag of features, called spatially sensitive bags of features (SS-BoF) [Bronstein:11]. The SS-BoF is a global descriptor defined in terms of mid-level features and the heat kernel, and can be represented by a square matrix whose elements represent the frequency of appearance of nearby codewords in the vocabulary. In the same spirit, Bu et al. [Bu:14] recently proposed the geodesic-aware bags of features (GA-BoF) for 3D shape classification by replacing the heat kernel in SS-BoF with a geodesic exponential kernel.
In this paper, we introduce a global spectral graph wavelet framework that characterizes the shape of the cortical surface of a carpal bone in terms of eigensystem of Laplace-Beltrami operator. In a bid to aggregate the local descriptors and transform them into global representation we multiply the global spectral graph wavelet matrix of each shape by its area matrix. The resultant global descriptor is not only isometric invariant, but also much more efficient and requires less memory storage. We prove through experiment on publicly-available database [Moore:07] that our proposed GSGW approach yields better performance in terms of MANOVA and permutation test compared to existing methods.
The rest of this paper is organized as follows. In Section 2, we briefly overview the Laplace-Beltrami operator. In Section 3, we introduce a global spectral graph wavelet framework, and we discuss in detail its main algorithmic steps. Experimental results are presented in Section 4. Finally, we conclude in Section 5 and point out some future work directions.
A 3D shape is usually modeled as a triangle mesh whose vertices are sampled from a Riemannian manifold. A triangle mesh may be defined as a graph or , where is the set of vertices, is the set of edges, and is the set of triangles. Each edge connects a pair of vertices . Two distinct vertices are adjacent (denoted by or simply ) if they are connected by an edge, i.e. .
2.1 Laplace-Beltrami Operator
Given a compact Riemannian manifold , the space of all smooth, square-integrable functions on is a Hilbert space endowed with inner product , for all , where (or simply ) denotes the measure from the area element of a Riemannian metric on . Given a twice-differentiable, real-valued function , the Laplace-Beltrami operator (LBO) is defined as , where is the intrinsic gradient vector field and is the divergence operator [Krim:15, Rosenberg:97]. The LBO is a linear, positive semi-definite operator acting on the space of real-valued functions defined on , and it is a generalization of the Laplace operator to non-Euclidean spaces.
Discretization A real-valued function defined on the mesh vertex set may be represented as an -dimensional vector , where the th component denotes the function value at the th vertex in . Using a mixed finite element/finite volume method on triangle meshes [Meyer:03], the value of at a vertex (or simply ) can be approximated using the cotangent weight scheme as follows:
where and are the angles and of two faces and that are adjacent to the edge , and is the area of the Voronoi cell at vertex . It should be noted that the cotangent weight scheme is numerically consistent and preserves several important properties of the continuous LBO, including symmetry and positive semi-definiteness [Wardetzky:07].
Spectral Analysis The matrix associated to the discrete approximation of the LBO is given by , where is a positive definite diagonal matrix (mass matrix), and is a sparse symmetric matrix (stiffness matrix). Each diagonal element is the area of the Voronoi cell at vertex , and the weights are given by
where and are the opposite angles of two triangles that are adjacent to the edge .
The eigenvalues and eigenvectors ofcan be found by solving the generalized eigenvalue problem using for instance the Arnoldi method of ARPACK, where are the eigenvalues and are the unknown associated eigenfunctions (i.e. eigenvectors which can be thought of as functions on the mesh vertices). We may sort the eigenvalues in ascending order as with associated orthonormal eigenfunctions , where the orthogonality of the eigenfunctions is defined in terms of the -inner product, i.e.
We may rewrite the generalized eigenvalue problem in matrix form as , where is an diagonal matrix with the on the diagonal, and is an orthogonal matrix whose th column is the unit-norm eigenvector .
The ability of eigenvectors of LBO for rendering the shape-based features is depicted in Figure 1. As shown, the lower-order eigenvectors capture the global structure of shape, while by increasing the number of eigenvectors more details of the curvature of the bone are captured.
The successful use of the LBO eigenvalues and eigenfunctions in shape analysis is largely attributed due to their isometry invariance and robustness to noise. Moreover, the eigenfunctions associated to the smallest eigenvalues capture well the large-scale properties of a shape. As shown in Figure 2, the (non-trivial) eigenfunctions of the LBO encode important information about the intrinsic global geometry of a shape. Notice that the eigenfunctions associated with larger eigenvalues oscillate more rapidly. Blue regions indicate negative values of the eigenfunctions and red colors regions indicate positive values, while green and yellow regions in between.
Figure 3 shows the normalized mean-squared error between the original bone surface and that reconstructed from an increasing number of eigenvectors of the LBO. As can be seen, just a small number of eigenvectors, i.e. , is sufficient to capture well features on the carpal bone surface to analyze shape differences in population study. Rendering carpal bone surface in lower-dimension enjoys some advantages like not being vulnerable to tessellation noise or image segmentation.
In this section, we provide a detailed description of our GSGW framework for analysis of cortical surface of a carpal bone using spectral graph wavelets. We start by defining the spectral graph wavelet transforms on a Riemannian manifold. We show how to build local descriptors from spectral graph wavelets and its subcomponent functions. Then, we propose a novel method to aggregate local descriptors to make global ones. Finally, we provide the main algorithmic steps of our carpal bone analysis framework.
3.1 Local Descriptors
Wavelets are useful in describing functions at different levels of resolution. To characterize the localized context around a mesh vertex , we assume that the signal on the mesh is a unit impulse function, that is at each mesh vertex . The spectral graph wavelet coefficients are expressed as
and that the coefficients of the scaling function are
Following the multiresolution analysis, the spectral graph wavelet and scaling function coefficients are collected to form the the spectral graph wavelet signature at vertex as follows:
where is the resolution parameter, and is the shape signature at resolution level given by
The wavelet scales () are selected to be logarithmically equispaced between maximum and minimum scales and , respectively. Thus, the resolution level determines the resolution of scales to modulate the spectrum. At resolution , the spectral graph wavelet signature is a 2-dimensional vector consisting of two elements: one element, , of spectral graph wavelet function coefficients and another element, , of scaling function coefficients. And at resolution , the spectral graph wavelet signature is a 5-dimensional vector consisting of five elements (four elements of spectral graph wavelet function coefficients and one element of scaling function coefficients). In general, the dimension of a spectral graph wavelet signature at vertex can be expressed in terms of the resolution as follows:
Hence, for a -dimensional signature , we define a spectral graph wavelet signature matrix as , where is the signature at vertex and is the number of mesh vertices. In our implementation, we used the Mexican hat wavelet as a kernel generating function . In addition, we used the scaling function given by
where and is set such that has the same value as the maximum value of . The maximum and minimum scales are set to and .
The geometry captured at each resolution of the spectral graph wavelet signature can be viewed as the area under the curve shown in Figure 4. For a given resolution , we can understand the information from a specific range of the spectrum as its associated areas under . As the resolution increases, the partition of spectrum becomes tighter, and thus a larger portion of the spectrum is highly weighted.
|(a) Mexican hat kernel for R=1||(b) Mexican hat kernel for R=5|
|(c) Mexican hat kernel for R=15||(d) Mexican hat kernel for R=30|
3.2 Global Descriptors
We refer to the -dimensional vector consisting of the first wavelet and scaling functions as the spectral graph wavelet point signature at vertex . Hence, we may represent the shape by an spectral graph wavelet signature matrix as of signatures, each of which is of length . In other words, the rows of are data points and the columns are features.
In a bid to aggregate the local descriptors and build global descriptor to represent a carpal bone surface we have to exploit bag-of-features paradigm. A major drawback of the BoF model is that it only considers the distribution of the codewords and disregards all information about the spatial relations between features, and hence the descriptive ability and discriminative power of the BoF paradigm may be negatively impacted. In addition, the BoF process is time-consuming since it requires different steps such as constructing dictionary, feature coding and feature pooling. To circumvent this limitation, we represent the shape by the vector , where is a area vector. We refer to as the global spectral graph wavelet (GSGW) descriptor of the carpal bone surface. In addition to circumvent the BoF process, the GSGW descriptor enjoys a number of desirable properties including simplicity, compactness, invariance to isometric deformations, and computational feasibility. Moreover, our spectral graph wavelet signature combines the advantages of both band-pass and low-pass filters.
3.3 Proposed Algorithm
Our proposed carpal bone analysis algorithm consists of two main steps. The first step is to represent each bone in the dataset by a spectral graph wavelet signature matrix, which is a feature matrix consisting of local descriptors. More specifically, let be a dataset of carpal bones modeled by triangle meshes . We represent each surface in the dataset by a spectral graph wavelet signature matrix , where is the -dimensional local descriptor at vertex and is the number of mesh vertices.
In the second step, the spectral graph wavelet signatures are aggregated to feature vector using computing for each carpal bone . Subsequently, all feature vectors of all shapes in the dataset are arranged into a data matrix . Figure 5 displays the GSGW codes of two carpal bones (capitate and lunate) from two different classes of wrist dataset [Moore:07] .
To assess the performance of the proposed GSGW framework, we employed two commonly-used evaluation criteria, the MANOVA and permutation test, which will be discussed in more detail in the next section.
4 Experimental Results
In this section, we evaluate the performance of our proposed GSGW approach for analysis of carpal bone surfaces through extensive experiments. The effectiveness of our method is validated by performing a comprehensive comparison with recent GPS approach [Chaudhari:14].
Datasets In order to evaluate the performance of the proposed GSGW framework on carpal bone surfaces a total of men and women with average age of years old from publicly-available benchmark [Moore:07] have been chosen (see Figure 6). More specifically, carpal bones of eight surfaces including capitate, hamate, lunate, pisiform, scaphoid, trapezium, trapezoid, triquetrum, and the first metacarpal for both wrists are used for performing analysis. Since the trapeziometacarpal joint of the thumb is a common site of osteoarthritis, the first metacarpal bone has been considered for our analysis. Each surface is rendered by triangular mesh with vertices and edges.
Performance Evaluation Measures In a bid to compare the shapes of the carpal bones in women versus men, we generated the GSGW signature for each carpal bone (eight total) and the first metacarpal bone for each subject for both the right and left wrists. For fair comparison of the performance of our proposed GSGW approach with GPS embedding, we followed setting in [Moore:07]
. we also compared the GSGW for each bone of the right and left wrist separately for the two groups (ten women versus ten men) using one-way multivariate analysis of variance (MANOVA) based on sex. Furthermore, in order to to further confirm the excellency of our approach we exploited non-parametric permutation tests as the resampling procedure.
For permutation testing, gender labels of the samples are randomly shuffled for
times, and test statistics calculated the to generate the null distribution. Then, we created the same statistics without shuffling and compute the percentile of that score, which referred to as the-value of the non-parametric permutation test. We report the -value acquired from both MANOVA and permutation testing. For , there would be a statistically significant difference between the two groups.
It should be noted that unlike GPS embedding in which the number of GPS coordinates varied for each carpal bone, we used a fixed dimension of which is acquired by computing principle component analysis of data matrix .
Baseline Method For the wrist benchmark [Moore:07] used for experimentation, we will report the comparison results of our method against GPS embedding approach [Chaudhari:14].
Implementation Details The experiments were conducted on a desktop computer with an Intel Core i5 processor running at 3.10 GHz and 8 GB RAM; and all the algorithms were implemented in MATLAB. The appropriate dimension (i.e. length or number of features) of a shape signature is problem-dependent and usually determined experimentally. For fair comparison, we used the same parameters that have been employed in the baseline methods, and in particular the dimensions of shape signatures. In our setup, a total of 31 eigenvalues and associated eigenfunctions of the LBO were computed. For the proposed approach, we set the resolution parameter to (i.e. the spectral graph wavelet signature matrix is of size , where is the number of mesh vertices).
4.1 Carpal Bone Dataset
Carpal bone dataset consists of mesh models from classes [Moore:07]. Each class contains objects with distinct postures. Moreover, each model in the dataset has approximately vertices.
Results In our GSGW approach, each surface in the carpal bone dataset is represented by a matrix of spectral graph wavelet signatures. The global spectral graph wavelet yields a vector of spectral graph wavelet codes, resulting in a data matrix of size . Figure 5 shows the spectral graph wavelet codes of two carpal bones (capitate and lunate) from two different classes of carpal bone dataset. As can be seen, the global descriptors are quite different and hence they may be used efficiently to discriminate between surfaces in statistical analysis tasks.
We compared our proposed GSGW method to GPS embedding method [Chaudhari:14] in terms of MANOVA and non-parametric permutation test. The results are summarized in Table 7 for right wrist and Table 8 for left wrist, respectively. Highlighted numbers in Tables show that the -value is exceeded .
As can be seen, our method achieves better analytical performance than GPS embedding method for both right and left wrist. For right wrist, the GSGW approach yields the lower -value compared with GPS embedding for six carpal bones out of nine ones in terms of MANOVA and for seven bones out of nine ones in terms of permutation test, respectively. In addition, the result of -value of MANOVA for some bones, e.g. hamate, has improved up to . For left wrist, our GSGW approach significantly improves the results by resulting the lower -value in terms of MANOVA for all bones except Metacarpal-1. Also, our method achieved much lower -value in terms of permutation test for all carpal bones in contrast with GPS embedding. Moreover, the result of -value of MANOVA for some bones, e.g. hamate, has improved up to . To speed-up experiments, all shape signatures were computed offline, albeit their computation is quite inexpensive due in large part to the fact that only a relatively small number of eigenvalues of the LBO need to be calculated.
In order to further assess the discriminative power of our approach, we computed the GSGW of carpal bone surfaces of the same class. As can be shown in Figures 9 and 10, even for very similar carpal bones with slightly difference, the proposed GSGW is able to distinguish and recognize the shape.
Parameter Sensitivity The proposed approach depends on two key parameters that affect its overall performance. The first parameter is the resolution of the spectral graph wavelet. The second parameter is the number of eigenvalues of LBO, which plays an important role in the GSGW vector . As shown in Figure 11, the best MANOVA result on carpal bones dataset is achieved using and . In addition, the performance of proposed method in terms of MANOVA is satisfactory for a wide range of parameter values, indicating the robustness of the proposed GSGW framework to the choice of these parameters.
In this paper, we presented a spectral-based framework for population study of carpal bones. Moreover, we performed statistical analysis on carpal bones of the human wrist by representing the cortical surface of the carpal bone using spectral graph wavelet descriptor to supply a means for comparing shapes of the carpal bones across populations. We also proposed a novel framework of directly extracting global descriptor so-called global spectral graph wavelet. Thus, we circumvented all the procedure of BoF paradigm which leads to a lower computation time as well as higher analysis accuracy. This approach not only captures the similarity between feature descriptors, but also substantially outperforms state-of-the-art methods both in accuracy and in scalability. For future work, we plan to apply the proposed approach to other surface analysis problems, and in particular segmentation and clustering.