The current study has two main findings. First, it is possible to analyze differences in surface morphometry by building a set of parametric surfaces using concepts from exterior calculus, such as differential one-forms and conformal nets. This is a high-level branch of mathematics that has not been extensively used in brain imaging before, but it provides a rigorous framework for representing, splitting, parameterizing, matching, and measuring surfaces. Second, the analysis of parametric meshes in computational studies of brain structure can be made more powerful by analyzing the full multivariate information on surface morphology. In this context, we are using multivariate to mean a multi-component description of surface shape at each and every point, rather than treating the observations at many different points on a surface as a multivariate vector (which can also be done, for example, in statistical shape analysis that use PCA on shapes). Work on tensor-based analysis of surfaces was proposed initially by Chung et al. (2003a)
, who noted that cortical surface data in children could be analyzed for patterns of growth over time by considering the local contraction or expansion of parametric grids adapted to the cortex (see Thompson et al. (1997
) for related work). A key barrier to getting these methods to work in general on subcortical surfaces is the complex topology of many subcortical surfaces, which this paper provides a method to understand.
In terms of validating the method, at least two types of validation are necessary. It is important to first show that the methods do in fact create conformal mappings (as checked by the assessment of angles in the resulting grids) and, second, show by formal proof that the patches form a total cover of the surface and are guaranteed not to overlap. For generating non-overlapping patches that are each proven to be conformal maps, the holomorphic one-form is a well-studied mathematical tool. According to Hodge theory, the holomorphic one-form induces a conformal mapping (i.e., parameterization) between a surface and the complex plane. The conformal parameterization induces a simple global structure, the conformal net. There are formal proofs that the resulting surface partition consists of non-overlapping patches that jointly cover the surface, and each patch is either a topological disk or a topological cylinder (Strebel, 1984
; Luo, 2006
). Basically, our segmentation algorithm just traces a critical graph so that each of the surface segments obtained is either a topological disk or a topological cylinder, and there is no overlap between patches.
In prior work, we also verified some of the formal properties of these maps, using artificially-generated, synthetic surfaces. In Wang et al. (2005b)
, we built a two-hole torus surface (essentially like a 3D solid version of a ), and found that the holomorphic one-form method was able to correctly split it into two patches with a rectangular parameterization. Via texture mapping, we also showed that the hippocampus can be represented as an open-boundary genus-one surface (a cylinder), and we built one-form based single-patch parameterizations. These correctly mapped the patch boundaries to the boundaries of the hippocampal surface, without causing the extremely dense gridding that is evident at the poles of a spherical parameterization. In Wang et al. (2005c)
, we used empirical plots to verify that the resulting maps were very close to conformal, even when internal landmark constraints were added to enforce landmark correspondences inside the conformal maps. Histograms of the angle difference from a conformal (angle-preserving) mapping were tightly clustered around zero even when a large number of landmark constraints were enforced (Wang et al., 2005c
; Lui et al., 2006a
). In Wang et al. (2005a)
, we also showed texture-mapped data verifying the surface registrations. In that paper, we used an additional term to improve the mutual information between hippocampal surface features (scalar fields defined on the surfaces) in two different subjects. Plots of the conformal factor and mean curvature were shown, before and after surface registration, to verify that the registration improved the alignment of corresponding features in parameter space, and, via the pull-back mapping, on the 3D surfaces. In the current work, we did not use mutual information to align curvature maps within the surface patches. This is because there is no reason to think that curvature is a reliable guide to homology on the ventricular surface instead, it seems reasonable to set up a computed correspondence between the two surfaces, that is smooth and conformal, and matches all the horns at once. If features on the ventricular surfaces were identified that could be used to independently validate the registration mappings, they could be used for verification, or included in the mapping themselves. However, there is no ground truth surface correspondence for ventricles across subjects, except for the logical requirement that each horn maps to its homolog (as is seen in our maps).
Figure 8 Detection of a synthetic group difference applied to the left ventricular surfaces of control subjects. The synthetic group of surfaces is generated by applying a geometric deformation (a known mathematical function) to a selected “bump” (more ...)
Finally, there may be rare pathological cases in which the patches in the critical graph may not correspond anatomically across different subjects, but we did not see any such examples in tests on large numbers of subjects. As the three horns on each side of the lateral ventricles are quite elongated, we found that in practice the patch boundaries always generated a cross on the surface (see ) in which the temporal and occipital horns met at a point, and the frontal horn wrapped around them in such a way that two points in the frontal horns parameter space met at the same 3D point. So long as the three tips of the ventricular horns can be isolated, such a decomposition is enforced, because the corners of the patches in parameter space contain right angles, so they still do when mapped into 3D via an angle-preserving (conformal) map.
We also assessed the sensitivity of the algorithm to simulated differences, by introducing spatially varying deformations of into the left ventricular surfaces of the control group. For each left ventricular surface, we selected a rectangular area in the parameter space of the frontal horn, as shown in . We slightly adjusted each surface point’s geometry coordinates by the following equation,
are tunable parameters that allowed us to vary the geometric deformation. Using surface registration, we introduced this geometric deformation to all the surfaces in the control group, to generate a new group of surfaces, matched to the original ones. We applied our multivariate tensor-based morphometry method to assess the geometric differences between these two groups. Using this scheme, abnormalities are deliberately created and mathematically defined, and the nature and context of the deformations can be systematically varied to determine the conditions that affect detection sensitivity.
illustrates the results of the experiments on these synthetic datasets. (a) illustrates the selected rectangle (labeled in yellow) on a specific left ventricular surface. We conducted 3 sets of experiments with different parameters: r = 0.005, k = 1 (), r = 0.005, k = 7 () and r = 0.01, k = 7 (). Although the imposed deformations were relatively small scale geometric perturbations involving only one coordinate, our algorithm still successfully picked up the areas that were significantly different by construction. For the frontal horn regions that involved the simulated deformation, the overall (corrected) significance of the group difference maps were p = 0.1642, 0.0248 and 0.0004, respectively, for the deformations of increasing severity.
These results illustrate the graded response of the multivariate tensor-based morphometry algorithm, in assessing deformations of spatially varying magnitude across the ventricular surfaces. This also demonstrates the effectiveness of our algorithm for picking up subtle morphometric deformations across surfaces.
In our experiment, an important step is the automatic partitioning of a complex 3D surface, i.e., we cut open a ventricular surface at its three extreme points. This turns the surface into a genus 0 surface with 2 open boundaries that is topologically equivalent to the surface in . Obviously, the original ventricular surface is a genus 0 surface. Theoretically it is topologically equivalent to a sphere. However, the concave shape, complex branching topology and extreme narrowness of the inferior and posterior horns make it extremely difficult to compute a meaningful regular mapping from a ventricular surface to a sphere (Friedel et al., 2005
compares our conformal parameterization results with spherical parameterization results from a robust genus-zero surface parameterization approach. Both of them are visualized by the texture mapping of a checker board. Our parameterization results are much more regular and uniform than the spherical parameterization results. The use of differential forms and cohomology theory helps to generate stable and accurate surface meshes, and a registration mapping that provides computed correspondences among different lateral ventricular surfaces.
Figure 9 Comparison of two lateral ventricular surface parameterization results. (a) is its conformal parameterization with the canonical holomorphic one-form; (b) is its spherical parameterization result (Friedel et al., 2005). Both parameterizations are visualized (more ...)
Our differential one-form method for modeling anatomical surfaces differs from the SPHARM (spherical harmonic) method (see, e.g., Styner et al. (2005
)) and the Laplace-Beltrami eigenfunction technique for analyzing surface shapes (Shi et al., 2009
). First, SPHARM or Laplace-Beltrami eigenfunction-based methods model geometric surfaces using a set of coefficients that weight the spherical harmonic basis functions or the Laplace-Beltrami eigenfunctions. These functions represent the entire shape of an anatomical structure as a sum of basis shapes (as we proposed in Thompson and Toga (1996
)). Instead, our method uses differential forms to generate global conformal surface parameterizations via mappings to a canonical parameter space, such as the complex plane here, for further surface registration and shape analysis. For our canonical conformal differential one-forms, all coefficients of the basis differential one-form functions are fixed and have a value of 1 (see Equation 2
). The surface is modeled by the one-form functions and not by their coefficients. Secondly, to combine or compare data across subjects, Styner et al. (2005)
) used parameter-space rotation based on the first-order term in the spherical harmonic expansion of the shapes, to register two anatomical model surfaces. Instead, we register surfaces by a constrained harmonic map in the canonical space, which is a smooth mapping with mathematically enforced regularity. Thirdly, we perform non-parametric permutation tests on tensors at each surface point, and assess the overall significance of group differences by non-parametric permutation tests on the area of statistics exceeding a pre-defined fixed threshold (Thompson et al., 2005b
). Since each tensor is defined based on local surface geometry, our map derives detailed local information on surface geometry, from the surface derivatives. In some cases, it may pick up more subtle differences than SPHARM, especially if the differences are highly localized and alter the local surface metric. In general, eigenfunction methods that generate shape coefficients (like SPHARM) will tend to be sensitive to more diffuse effects on the global shape of a structure, especially ones that can be represented accurately by a projection onto the first few eigenfunctions of the associated surface differential operator, whereas our map-based analysis will also tend to be sensitive to highly localized differences in surface geometry. Eigenfunction methods may miss small, localized effects, unless a large number of basis functions is used. In addition, it may be more intuitive in brain mapping to show statistical differences as a 3D map, rather than as a difference in shape coefficients, which may not have any convenient anatomical interpretation.
In the most related work, Chung et al. (2003a
have also proposed to use the surface metrics as the basis for morphometry in cortical studies. From the Riemannian metric tensors, they computed the local area element. They further defined the surface area dilatation, which is approximately the trace of the Jacobian determinant of the mapping from the template to the cortical surface. At each point on the surface, they used a Student’s t
test to assess group differences between patients with autism and matched control subjects. Although we used a different cortical surface registration algorithm, their proposed statistic is in fact very close to the determinant of Jacobian matrix, which we evaluated in our comparison experiments (). The new method we propose in this paper detects group difference in the “Log-Euclidean” space, i.e., using Lie group methods to handle the curvature of the tensor manifold of symmetric positive definite matrices. As illustrated in our experimental comparisons, the new method retains full Riemannian metric information and may capture additional useful information on surface morphometry.
The work reported here is related to ongoing research on the conformal mapping of brain surfaces. Brain surface conformal parameterization has been studied intensively (Schwartz et al., 1989
; Balasubramanian et al., 2009
; Hurdal and Stephenson, 2004
; Angenent et al., 1999
; Gu et al., 2004
; Ju et al., 2005
; Wang et al., 2006
). Most brain conformal parameterization methods (Hurdal and Stephenson, 2004
; Angenent et al., 1999
; Gu et al., 2004
; Ju et al., 2004
; Joshi et al., 2004
; Ju et al., 2005
) can handle the entire cortical surface of the brain, but cannot usually handle surfaces with boundaries or junctions with other surfaces. The canonical holomorphic one-form method (Wang et al., 2007
), the Ricci flow method (Wang et al., 2006
) and slit map method (Wang et al., 2008a
) are ideal for these situations, as they can handle surfaces with complicated topologies.
The multivariate TBM method proposed here for analyzing surfaces is quite general. In most engineering studies, surfaces are usually represented by triangulated meshes. Once the surface mapping is established, the Jacobian matrix can always be computed from Equation 3
, and multivariate TBM can then be directly applied. The only difference between this and conventional TBM is that logarithmic transforms are applied to convert the tensors into vectors that are more tractable for use with Euclidean operations, and Hotelling’s T2
test is applied because the resulting data is multivariate at each point. Besides the proposed constrained harmonic map, other popular surface registration methods such as CARET (Van Essen et al., 2001
), BrainVoyager QX (http://www.BrainVoyager.com
), FreeSurfer (Fischl et al., 1999
), and the cortical pattern matching method (Thompson et al., 2004b
) can also generate parametric surfaces as inputs for our multivariate TBM method.
When tensor-based morphometry (TBM) is applied to 3D images, the interpretation of the determinant of the Jacobian is obvious: it represents an estimate of the volume difference between two images, such as a patient image and a template, or between two successive images of the same subject. In our analysis, we also derived other metrics from the Jacobian (or, strictly speaking, from the surface metric tensor). In particular, the eigenvalues of the surface metric tensor reflect a tensors dilation effects in two orthogonal directions. They carry partial information about the tensor, and may pick up on changes that occur along one direction in the surface grid. For instance, in Leporé et al. (2008)
and Brun et al. (2009a)
, we argued that, in general, if a structure is stretched preferentially in one direction, but contracted in another, a tensor-based analysis is likely to pick up this effect, but a volumetric analysis may miss it. A more plausible case is that a disease that causes atrophy in all directions, but slightly greater atrophy in one direction, In Alzheimer’s disease, for example, the cross-sectional area of the hippocampus may shrink more rapidly than the anterior-posterior length of the hippocampus. In the current analysis, if the dilation of the ventricles were more pronounced in a radial than an anterior-posterior direction, or, if the normal variations in the ventricles were greater in a radial than an anterior-posterior direction, then the eigenvalues of the surface tensor could be used to pick up these differential effects with high sensitivity. For instance, the biological meaning of a statistical difference in a tensor eigenvalue would be that the structure was abnormally dilated (on contracted) along the grid lines on the surface mesh, which may reflect a radial expansion, or an anterior-posterior stretching, or both. We admit that in general, a difference in the tensors may initially be more difficult to interpret than a difference in a structures volume or surface area. Even so, multivariate TBM detecting a group difference could be followed up with an analysis of the individual eigenvalues or rotations in the tensor to hone in on a more specific biological interpretation of the salient differences.
As we noted in Leporé et al. (2008)
, there are two main caveats when applying this multivariate version of TBM. First, just because a method detects group differences with greater effect sizes, it does not mean that they are true, as we do not have ground truth information on the true extent of ventricular expansion in HIV. Even so, it is logical from observing the maps in that any such group differences in surface tensors arise from true systematic differences in brain structure that lead to distortions in surface morphology. Second, it is not always the case that the multivariate version of TBM must give higher effect sizes than the standard univariate version of TBM (i.e., analysis of group differences in regional volumes encoded in the Jacobian determinant). We have found in the past that when sample sizes are relatively small, the covariance terms that incorporate empirical information on brain variation have a large number of free parameters, needing relatively large numbers of subjects estimate them accurately (Brun et al., 2009b
Even so, we were able to show the power advantage of the multivariate version of TBM even in this very small sample of subjects with ventricular surface data. This is partly because the differences in mean ventricular volume are extremely large in HIV patients versus controls (Thompson et al., 2006
), but it shows the potential of the approach for use in small samples.
In future, we will apply our multivariate TBM framework to additional 3D MRI datasets to study brain surface morphometry. We plan to apply other holomorphic differentials, such as holomorphic quadratic differentials, and meromorphic differential forms to study related surface regularization, parameterization and registration problems.