Statistical Shape Analysis of Brain Structures using SPHARM-PDM

— Shape analysis has become of increasing interest to the neuroimaging community due to its potential to precisely locate morphological changes between healthy and pathological structures. This manuscript presents a comprehensive set of tools for the computation of 3D structural statistical shape analysis. It has been applied in several studies on brain morphometry, but can potentially be employed in other 3D shape problems. Its main limitations is the necessity of spherical topology. The input of the proposed shape analysis is a set of binary segmentation of a single brain structure, such as the hippocampus or caudate. These segmentations are converted into a corresponding spherical harmonic description (SPHARM), which is then sampled into a triangulated surfaces (SPHARM-PDM). After alignment, diﬀerences be-tween groups of surfaces are computed using the Hotelling T 2 two sample metric. Statistical p-values, both raw and corrected for multiple comparisons, result in signiﬁcance maps. Additional visualization of the group tests are provided via mean diﬀerence magnitude and vector maps, as well as maps of the group covariance information. The correction for multiple comparisons is performed via two separate methods that each have a distinct view of the problem. The ﬁrst one aims to control the family-wise error rate (FWER) or false-positives via the extrema histogram of non-parametric permutations. The second method controls the false discovery rate and results in a less conservative estimate of the false-negatives.


I. Introduction
Quantitative morphologic assessment of individual brain structures is often based on volumetric measurements.Volume changes are intuitive features as they might explain atrophy or dilation due to illness.On the other hand, structural changes at specific locations are not sufficiently reflected in volume measurements.Shape analysis has thus become of increasing interest to the neuroimaging community due to its potential to precisely locate morphological changes between healthy and pathological structures.
One of the first and most influential research in shape analysis was presented by D'Arcy [1] in his ground-breaking book On Growth and Form.In more recent years, several researchers proposed shape analysis via deformable registration to a template ( [2], [3], [4], [5]).Inter-subject comparisons are made by analyzing the individual deformable transformations.This analysis of the transformation fields has to cope with the high dimensionality of the transformation, the template selection problem and the sensitivity We are thankful to Steve Pizer and Sarang Joshi for highly valuable discussions in regard to shape analysis.This research has been/is supported by the UNC Neurodevelopmental Disorders Research Center HD 03110, NIH Roadmap for Medical Research, Grant U54 EB005149-01, NIH NIBIB grant P01 EB002779, and the EC-funded BIOMORPH project 95-0845.
to the initial position.Nevertheless, several studies have shown stable shape analysis results.[6] and [7] presented some of the first mathematical methods for 3D shape analysis based on sampled descriptions.The shape analysis of densely sampled 3D Point Distribution Models (PDM) and their deformations was first investigated by [8].Inspired by their experiments, [9] proposed shape analysis based on a parametric boundary description called SPHARM ( [10]).The SPHARM shape analysis approach was extended by [11] to use the implied sampled surface (SPHARM-PDM), a method also used by [12].[13], [14] and Golland [15] proposed shape analysis on medial shape descriptions in 3D and 2D, respectively.They used a fixed topology sampled model with implicit correspondence that is fitted to the objects.
This manuscript presents a comprehensive set of tools that form a pipeline for the computation of 3D structural statistical shape analysis on the object boundary via SPHARM (see Figure 1).It has been applied in several studies on brain morphometry, but can potentially be employed in other 3D shape analysis problems.The main limitations of our shape analysis methodology is the necessity of spherical topology.

II. Methods
This section describes the SPHARM-PDM shape description, the statistical analysis methodology, as well as issue of alignment and scaling normalization.
In summary, the input of the proposed shape analysis is a set of binary segmentation of a single brain structure, such as the hippocampus or caudate.These segmentations are first processed to fill any interior holes and a minimal smoothing operation.The processed binary segmentations are converted to surface meshes, and a spherical parametrization is computed for the surface meshes using a area-preserving, distortion minimizing spherical mapping.The SPHARM description is computed from the mesh and its spherical parametrization.Using the first order ellipsoid from the spherical harmonic coefficients, the spherical parametrizations are aligned to establish correspondence across all surfaces.The SPHARM description is then sampled into a triangulated surfaces (SPHARM-PDM) via icosahedron subdivision of the spherical parametrization.These SPHARM-PDM surfaces are all spatially aligned using rigid Procrustes alignment.Group differences between groups of surfaces are computed using the standard robust Hotelling T 2 two sample metric.Statistical p-values, both raw and corrected for multiple comparisons, result in significance maps.Additional visualization of the group tests are provided via mean difference magnitude and vector maps, as well as maps of the group covariance information.
The correction for multiple comparisons is performed via two separate methods that each have a distinct view of the problem.The first one aims to control the family-wise error rate (FWER) or false-positives via the extrema histogram of non-parametric permutations.The second method controls the false discovery rate and results in a less conservative estimate of the false-negatives.In summary, the SPHARM description is a hierarchical, global, multi-scale boundary description that can only represent objects of spherical topology ( [10]).The spherical parameterization is computed via optimizing an equal area mapping of the 3D voxel mesh onto the sphere and minimizing angular distortions.The basis functions of the parameterized surface are spherical harmonics.Each individual SPHARM description is composed of a set of coefficients, weighting the basis functions.[16] demonstrated that SPHARM can be used to express shape deformations.Truncating the spherical harmonic series at different degrees results in object representations at different lev-els of detail.SPHARM is a smooth, accurate fine-scale shape representation, given a sufficiently high representation level.Based on a uniform icosahedron-subdivision of the spherical parameterization, we obtain a Point Distribution Model (PDM).

A.1 Mathematics behind SPHARM
This section discusses the mathematical properties of the spherical harmonic basis functions.It gives a summary of spherical harmonic descriptors as they are presented in Brechbühler's dissertation [17].
Spherical harmonic basis functions Y m l , −l ≤ m ≤ l of degree l and order m are defined on θ ∈ [ 0; π ] × φ ∈ [ 0; 2π) by the following definitions: where Y m l * denotes the complex conjugate of Y m l and P m l the associated Legendre polynomials Table I lists explicit expressions for the spherical harmonic functions up to degree 3. The Cartesian notion reveals that the spherical harmonics are polynomials in the 3D space (u 0 , u 1 , u 2 ).To express a surface using spherical harmonics, the three coordinate functions are decomposed and the surface v(θ, φ) = (x(θ, φ), y(θ, φ), z(θ, φ)) T takes the form where the coefficients c m l are three-dimensional vectors due to the three coordinate functions.The coefficients c m l are obtained by solving a least-squares problem.Therefore, the values of the basis functions are gathered in the matrix Z = (z i,j(l,m) ) with z i,j(l,m) = Y m l (θ i , φ i ), where j(l, m) is a function assigning an index to every pair (l, m) and i denotes the indices of the n vert points to be approximated.The coordinates of these points are arranged in V = ( v 1 , v 2 , . . ., v nvert ) T and all coefficients are gathered in C = ( c 0 0 , c −1 1 , c 0 1 , . ..)T .The coefficients that best approximate the points in a least-squares sense are obtained by Using spherical harmonic basis functions, we obtain a hierarchical surface description that includes further details as more coefficients are considered.This is illustrated in Fig. 2.

A.2 SPHARM Correspondence
Correspondence of SPHARM-PDM is determined by normalizing the alignment of the spherical parameterization to an object-specific frame.In the studies presented in this paper, the normalization is achieved by rotation of the parameterization, such that the spherical equator, 0 • and 90 • longitudes coincide with those of the first order ellipsoid(see Figure 3).After this parametrization normalization, corresponding surface points across different objects possess the same parameterization.The SPHARM-PDM shape analysis is visualized in Figure 1 using a caudate structure.Prior to the shape analysis, the group average object is computed for each subject group, and an overall average object is computed over all group average objects.Each average structure is computed by averaging the 3D coordinates of corresponding surface points across the group.The overall average object is then used in the shape analysis as the template object.At every boundary point for each object, we compute a distance map representing the signed local Euclidean surface distance to the template object.The sign of the local distance is computed using the direction of the template surface normal.In the global shape analysis, the average of the local distances across the whole surface is analyzed with a standard group mean difference test.The local shape analysis is computed by testing the local distances at every boundary point.This results in a significance map that represents the significance of these local statistical tests and thus allows locating significant shape differences between the groups.We corrected the shape analysis for the multiple comparison problem using a uniformly sensitive, non-parametric permutation test approach ( [18]).The non-corrected significance map is an optimistic estimate of the real significance, whereas the corrected significance map is a pessimistic estimate that is guaranteed to control the rate of false positives at the given levelα (commonly α = 0.05) across the whole surface.
We have compared the SPHARM correspondence to other types of correspondences and it compared favourably to human expert established correspondence [19].

B. Area-preserving Spherical Mapping
The appropriate parameterization of the points of a surface description is a key problem for correspondence finding, as well as for an efficient SPHARM representation.Every point i on the surface is to be assigned a parameter vector (θ i , φ i ) that are located on the unit sphere.A homogeneous distribution of the parameter space is essential for an efficient decomposition of the surface into SPHARM coefficients.We give here a brief summary of the surface parameterization procedure proposed by Brechbühler [10].
A bijective mapping of the surface to the unit sphere is created, i. e., every point on the surface is mapped to exactly one point on the sphere, and vice versa.The main idea of the procedure is to start with an initial parameterization.This initial parameterization is optimized so that every surface patch gets assigned an area in parameter space that is proportional to its area in object space.The initial mapping from surface to parameter space is constructed using discrete Laplace's equations to solve the corresponding Dirichlet problem.To obtain a homogeneous distribution of the parameter space over the surface, this initial parameterization is modified in a constrained optimization procedure considering two criteria: 1. Area preservation: Every object region must map to a region of proportional area in parameter space.2. Minimal distortion: Every quadrilateral should map to a spherical quadrilateral in parameter space.
Brechbühler establishes constraints for area preservation, while the distortion of the mesh serves as the objective function during optimization.The optimization solves the resulting system of nonlinear equations by linearizing them and taking Newton steps.

C. Surface Models from SPHARM: SPHARM-PDM
From the SPHARM description we can compute triangulated surfaces by sampling the spherical parameterization uniformly on the sphere.Equidistant sampling in the parameter space leads to a dense sampling around the poles (θ = 0, θ = π) and a coarse sampling around the equator (θ = π/2).This fact can be explained by the poles being mapped to all points having φ = 0 . . .2π and θ = 0, or θ = π (see also Fig. 3).Using a uniform icosahedron subdivision shown in Fig. 4, however, we gain a good approximation of a homogeneous sampling of the spherical parameter space and thus also of the object space.Using the pre-computed parameter locations (θ i , φ i ) from the icosahedron subdivision, we can compute the PDM directly from the coefficients, as the parameter locations stay constant for all objects at a given subdivision level.The sampling points x i at the locations (θ i , φ i ) are obtained using equation ( 4): In our experience the sampling for the surface of brain structures should be chosen between subdivision level 10 (1002 surface points) and 20 (4002 surface points) depending on the complexity of the objects.For example, hippocampal surfaces are well represented by a subdivision level 10 and caudate surfaces by a subdivision level 15.

D. Alignment and Scaling Normalization
Alignment and scaling of the objects are two important issues in shape analysis that are also discussed here [11].In general, we usually chose rigid-body Procrustes alignment ( [20]) followed by no scaling normalization or scaling inversely to the intra-cranial cavity volume (ICV).Commonly both results are reported in a shape analysis study, the original scale analysis represents an unbiased raw analysis, whereas the ICV scale analysis corrects for differences in overall interior head size and thus approximatively normalizes for gender and age differences.In general both types of analysis should yield similar patterns of differences and special care needs to be taken if this is not the case.

E. Local Testing of Group Mean Difference using Hotelling
T 2 metric After shape description, correspondence establishment, alignment and scaling normalization, the next step in the shape analysis pipeline consists of testing for differences between groups at every surface location (see Figure 5).This can be done in 2 main fashions: • Analyzing the magnitude of the local difference vector to a template: For this option, a template needs to be first selected, usually this is the common mean of the 2 groups or the mean of separate control group.The magnitude of the difference is easily understood and results in difference maps for each subject.Also, the resulting univariate statistical analysis is quite well known and local significance can be easily computed using the Student t metric.The main disadvantage of this method is the need to select a template, which introduces an additional bias into the statistical analysis.We applied this method in earlier studies of hippocampi [21] and ventricles [22].
• Analyzing the spatial location of each point: For this option, no template is necessary and multivariate statistics of the (x, y, z) location is necessary.We have chosen to use the Hotelling T 2 two sample difference metric as a measurement of how 2 groups locally differ from each other.The standard Hotelling T 2 is defined as is the pooled covariance matrix.An alternative modified Hotelling T 2 metric is less sensitive to group differences of the covariance matrixes and the number of samples [23]: . All our current studies are based on this modified Hotelling T 2 metric.
We then want to test the two groups for differences in the means of the selected difference metric (univariate: Student t, multivariate: Hotelling T 2 ) at each spatial location.Permutation tests are a valid and tractable approach for such an application, as they rely on minimal assumptions and can be applied even when the assumptions of the parametric approach are untenable.Non-parametric permutation tests are exact, distribution free and adaptive to underlying correlation patterns in the data.Further, they are conceptually straightforward and, with recent improvements in computing power, are computationally tractable.
Our null hypothesis is that the distribution of the locations at each spatial element is the same for every subject regardless of the group.Permutations among the two groups satisfy the exchangeability condition, i.e. they leave the distribution of the statistic of interest unaltered under the null hypothesis.Given n 1 members of the first group a k , k = 1 . . .n 1 and n 2 members of the second group b k , k = 1 . . .n 2 , we can create M ≤ n1 + n2 n2 permutation samples.A value of M from 20000 and up should yield results that are negligibly different from using all permutations [24].

F. Correction for Multiple Comparison
The local shape analysis involves testing from a few to many thousands of hypothesis (one per surface element) for statistically significant effects.It is thus important to control for the multiple testing problem (see Figure 5), and the most common measure of multiple false positives is the familywise error rate (FWER).The multiple testing problem has been an active area of research in the functional neuroimaging community.The most widely used methods in the analysis of neuroimaging data use random field theory [25] [26] and make inferences based on the maximum distribution.In this framework, a closed form approximation for the tail of the maximum distribution is available, based on the expected value of the Euler characteristic of the thresholded image [26].However, random field theory relies on many assumptions such as the same parametric distribution at each spatial location, a point spread function with two derivatives at the origin, sufficient smoothness to justify the application of the continuous random field theory, and a sufficiently high threshold for the asymptotic results to be accurate.
In our work, we are employing non-parametric permutation tests and false discovery rate as two alternative correction methods for the multiple comparison problem.

F.1 Minimum over Non-parametric permutation tests
Non-parametric methods deal with the multiple comparisons problem [18] and can be applied when the assumptions of the parametric approach are untenable.Similar methods have been applied in a wide range of functional imaging applications [27], [28], [29].
Details of this method are described in [18] and only a summary is provided here.The correction method is based on computing first the local p-values using permutation tests.The minimum of these p-values across the surface is then computed for every permutation.The appropriate corrected p-value at level α can then be obtained by the computing the value at the α-quantile in the histogram of these minimum values.Using the minimum statistic of the p-values, this method correctly controls for the FWER, or the false positives, but no control of the false negatives is provided.The resulting corrected local significance values can thus be regarded as pessimistic estimates akin to a simple Bonferroni correctection.

G. False Discovery Rate
Additionally to the non-parametric permutation correction, we have also implemented and applied a False Discovery Rate Estimation (FDR) [30], [31] method.The innovation of this procedure is that it controls the expected proportion of false positives only among those tests for which a local significance has been detected.The FDR method thus allows an expected proportion (usually 5%) of the FDR corrected significance values to be falsely positive.The correction using FDR provides an interpretable and adaptive criterion with higher power than the non-parametric permutation tests.FDR is further simple to implement and computationally efficient even for large datasets.
The FDR correction is computed as follows [31]: 1. Select the desired FDR bound q, e.g.5%.This is the maximum proportion of false positives among the significant tests that you are willing to tolerate (on average).2. Sort the p-values smallest to largest.3. Let p q be the p-value for the largest index i of the sorted p-values p sort,i ≤ q •i/N , where N is the number of vertices.4. Declare all locations with a p-value p ≤ p q significant.

III. Computation scheme and Implementation
This section describes the set of shape analysis tools and the use of these in succession to compute the whole shape analysis.
General overview with schematic view of how to use the tools (input/output) Detailed description of tools

IV. Results
Example results on the male caudate data

Fig. 3 .
Fig. 3. Left: Illustration of SPHARM correspondence via alignment of the spherical parametrization using the first order ellipsoid meridian and equator.Right: Set of 6 caudate structures with correspondence shown at selected locations via colored spheres.

Fig. 5 .
Fig. 5. Left: Schematic view of testing procedure: A set of features per surface point (shown as colormaps) are analyzed using unior multi-variate statistics, resulting in a raw significance map (bottom) .Right: The raw significance maps are overly optimistic and need to be corrected to control for the multiple comparison problem.The corrected significance maps on the other hand are commonly pessimistic estimates.

Fig. 6 .
Fig. 6.Top: Scheme of P-value computation via permutation tests.The real group difference S 0 (e.g.Hotelling T 2 ) is compared to that the group differences S j of random permutations of the group labels.The quantile in the S j histogram associated with S 0 is the p-value.Bottom: Example run on test data with 1000 random permutations.On the left the individual S j are plotted and on the left the corresponding S j histogram is shown.The black line indicates the value of S 0 .

TABLE I Explicit
expressions of the spherical harmonics up to degree 3, in both polar and Cartesian form due to Brechbühler.The last part of the table gives the common normalizing constants, e.g.Y 0 1 = 3/4πu 2 .