PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Proc IEEE Int Symp Biomed Imaging. Author manuscript; available in PMC 2010 August 17.
Published in final edited form as:
Proc IEEE Int Symp Biomed Imaging. 2004 April 1; 1: 428–431.
doi:  10.1109/ISBI.2004.1398566
PMCID: PMC2922987
NIHMSID: NIHMS224685

CORTICAL SURFACE PARAMETERIZATION BY P-HARMONIC ENERGY MINIMIZATION

Abstract

Cortical surface parameterization has several applications in visualization and analysis of the brain surface. Here we propose a scheme for parameterizing the surface of the cerebral cortex. The parameterization is formulated as the minimization of an energy functional in the pth norm. A numerical method for obtaining the solution is also presented. Brain surfaces from multiple subjects are brought into common parameter space using the scheme. 3D spatial averages of the cortical surfaces are generated by using the correspondences induced by common parameter space.

1. INTRODUCTION

The surface area of the gray matter of the cerebral cortex is approximately 1570 cm2 [8]. 60–70% of the surface area is buried in sulcul folds and creases. There is considerable variability in the size, location and extent of the sulci and gyri across human subjects. Bringing multiple brain surfaces into a common coordinate system is helpful in studying variability of these sulcal patterns across subjects, for integrating and averaging functional data across subjects, and in studying patterns in cortical development over time. Since the cortex can be modeled as a convoluted sheet with spherical topology, it is natural to parameterize it using spherical coordinates [7, 9].

Eck et al. [3] and Kanai et al. [4] model a triangulated surface as a configuration of springs with one spring placed along each edge of each triangle. The resulting energy functional, the harmonic energy, is shown to be a quadratic form and is minimized using gradient descent to transform the surface into a planar disk.

Desbrun et al. [5] propose a parameterization technique which uses cotangent of angles in the given triangulation. The resulting energy functional (Tuette energy) is argued to be a measure of angle distortion and its minimization leads to a new parameterization.

Haker et al. [2] presented a method for conformally mapping the cortical surface to a sphere. Their method uses the Laplace-Beltrami operator and the fact that for a conformal map, the Laplace-Beltrami operator calculated on the parameterization function is zero everywhere on the surface. Though these methods ensure a perfectly conformal map, the stereographic projections involved can introduce a large amount of length and area distortion.

Circle packing was introduced as a quasi-conformal parameterization method in [1]. Analytic surfaces can be approximated by circle packing, but the surface packing method considers only the connectivity and not the geometry of the surface [12, 10].

In this paper, a parameterization technique for the cortical surface based on p-harmonic energy minimization is discussed. We argue that since higher order norms result in smoother solutions, the resulting maps are smoother and more uniform. The problem is formulated and solved using Finite Elements.

2. MATHEMATICAL FORMULATION

The parameterization of an arbitrary surface, tesselated using a set of triangles, can be viewed as an assignment of complex numbers or vectors in R2 to each vertex in the triangulation. We assign a complex number to every point in the triangulation in such a way that the resulting p-harmonic energy is minimized. Let S be a submanifold of Rn of dimension 2. We define [var phi]: SR2 to be a function such that the p-harmonic energy functional is minimized.

The p-harmonic energy functional is defined as

Es=||φ||pdS
(1)

Here [var phi] is the unknown vector valued function defined over the given surface. The operator [nabla] represents the covariant derivative which can be viewed as the gradient with respect to the local coordinate system as discussed later.

We can see that if no constraints are imposed, then the minimization will result in [var phi] = 0 everywhere. Thus, appropriate constraints need to be imposed to find a non-trivial solution. Here we constrain the mapping of the inter-hemispheric fissure.

For p = 2 the energy functional Es is equivalent to the Dirichlet energy described in the conformal mapping literature. However if we impose more than two point constraints on the minimization, the resulting map will not be conformal. Eells-Sampson [13] prove the existence of p-harmonic maps in the case when the target manifold is a subset of Rn that has non-positive sectional curvature. The disk R2 in has these properties.

We write the energy functional as the sum of two energy functionals such that the corresponding arguments are scalars,

Es=||α||pdS+||β||pdS

where [var phi] = [α, β]′. Since α and β can be chosen independently, we can write,

minφEs=minφ||φ||pdS=minα||α||pdS+minβ||β||pdS

Therefore, this minimization can be performed by minimizing over the two real valued functions separately. Here we demonstrate minimization of one real-valued function using finite elements; minimization of the other is done in exactly the same way other than the constraints applied when mapping the inter-hemispheric fissure.

3. FINITE ELEMENT FORMULATION

The above real-valued functions are defined on the continuous domain which we discretize using finite elements. We assume that both α and β are piecewise linear functions defined on a triangular tessellation of the cortical surface.

To discretize the problem, we use a local coordinate system for each triangle. The use of a local coordinate system for each triangle is justified because we are only interested in the magnitude of the gradient and not its direction. The z-axis is chosen normal to the triangle and the x-axis lies on the longest edge of the triangle.

Let α be a piecewise linear real-valued scalar function defined over the surface, and αi is the function restricted to triangle i. Since αi is linear on the ith triangle we can write,

αi(x,y)=a0i+a1ix+a2iy
(2)

The three coefficients can be determined if values of the function α are known at the three vertices of the triangle. These equations can be written in matrix form as,

[1x1iy1i1x2iy2i1x3iy3i]Di[a0ia1ia2i]=[αi(x1,y1)αi(x2,y2)αi(x3,y3)]
(3)

The coefficients a0i,a1i and a2i can be obtained by inverting the matrix. From eq. 2 and by inverting the matrix in eq. 3, we obtain

αi=[a1ia2i]=1Di[y2iy1iy3iy1iy1iy2ix1ix2ix1ix3ix2ix1i]Bi[αi(x1,y1)αi(x2,y2)αi(x3,y3)]Γiαi=12AiBiΓi

We use the fact that for any triangle |Di| = 2Ai, where Ai is the area of the triangle. Since αi is piecewise linear, its gradient is constant over each triangle i, so that:

||α||pdS=i||αi||pAi

where the sum is over all triangles.

minα||α||pdS=minΓi(Ai)1p12BiΓip=minΓi||MiΓi||p=minΓ||MΓ||p

where Mi=(Ai)1p12Bi, M is composed using coefficients of Mi, and Γ is a vector with coefficients α for each vertex. The vector MΓ can be split into two parts: free vertices and constrained vertices. Values of α at constrained vertices are known.

minα||α||pdS=minΓf||MfΓf+McΓc||p

where Mf, Γf and Mc, Γc, are the free and constrained parts, respectively of the M and Γ matrices.

By removing the constrained vertices, we now have an unconstrained minimization problem. For p = 2, the solution is given by taking the pseudoinverse of the matrix Mf.

Γf=MfMcΓc
(4)

The fact that matrix Mf is sparse allows us to use the computationally efficient conjugate gradient method for obtaining the solution. The Jacobi preconditioner reduces the execution time considerably. For p > 2, we use a nonlinear conjugate gradient method. As noted above, we consider only even values of p so that the gradients are analytic. This is particularly easy for even values of p.

4. PLANAR AND SPHERICAL MAPPINGS

The goal of this work is to develop a common parameterization across multiple brains. While the parameterization applies to the original brain surface, for the purposes of illustration and visualization it is often convenient to represent these surfaces as mappings to a plane or sphere. To compute the parameterizations we first generated cortical surfaces using our BrainSuite software package [11] (Shattuck and Leahy, 2002), which includes a six-stage cortex modeling sequence. First, the brain is extracted from surrounding skull and scalp tissues using a combination of edge detection and mathematical morphology. Next, the intensities of the MRI are corrected for shading artifacts. Each voxel in the corrected image is labeled according to tissue type using a statistical classifier. The white matter in the brain volume is selected to produce a binary mask; the user then cuts through the brainstem in this mask to remove the brainstem and cerebellum. The resulting mask represents the cerebrum, but it is likely that tessellation of this volume will produce surfaces with topological handles. Prior to tessellation, these handles are identified and removed using a graph-based approach. The resulting mask is then tessellated to produce a genus zero surface.

Our spherical mapping is computed from the genus-zero surface in two stages. First, each brain hemisphere is mapped onto a unit disk. This is done by constraining a closed contour, representing the intersection of a plane through the inter-hemispheric fissure with the corpus callosum, to lie on the unit circle. We parameterize this contour on the unit circle by arc length. With this constraint we solve for the p-harmonic map, as described above, for each hemisphere in turn. We then use a stereographic projection to map the unit disk onto a hemisphere. Performing this operation for both hemispheres produces a bijective mapping of the cortical surface to the unit sphere.

5. EXAMPLES

An example of a single cortical surface and its spherical maps for p = 2, 4 and 6 are shown in Fig. 1. Note that as p increases, the mapping of cortical features becomes more uniform. A desirable property of the paramaterization is that the associated spherical or planar maps are minimally distorted and that corresponding gyral and sulcal features map to similar coordinates across subjects. This can significantly reduce the difficulty of the subsequent problem of automatically matching or identifying cortical features across a population.

Fig. 1
A cortical surface extracted from MRI data and its maps on a sphere for different values of p

As a preliminary evaluation of this method, we have investigated the average behavior of these maps for several subjects as a function of p. Volumetric MR images of six volunteers were co-registered using a 12-parameter affine fit with AIR [14]. Cortices were extracted and parameterized as described above. The maps were resampled onto a regular lattice with respect to the parameter [var phi] and partially aligned by rotation about the unit circle corresponding to the inter-hemispheric fissure. Shown in Fig. 2 are brain surfaces produced by averaging the original 3D coordinates for each parameterized brain at corresponding points on the regular lattice in [var phi].

Fig. 2
The figure shows spatial averages calculated for different values of p

Ideally the parameterizations would align sulcal landmarks across subjects so that the averaged brain would serve as an atlas containing clearly identifiable landmarks. Since the parameterization uses only knowledge of the location of the inter-hemispheric fissure, we do not expect to achieve this degree of alignment. However, the results in Fig. 2 appear to show that as p increases, and distortion decreases, there is indeed improved correspondence between brains so that the average reveals a more realistic representation of common cortical features.

We plan a more detailed evaluation of this approach in which we will investigate the impact of the value of p, and of the choice of volumetric alignment method prior to cortex extraction, on the degree of metric distortion required to bring cortices into alignment using covariant flows constrained by hand-labeled sulcal features [15].

References

1. Hurdal MK, Stephenson K, Bowers PL, Sumners DWL, Rottenberg DA. Coordinate systems for conformal cerebellar flat maps. NeuroImage. 2000;11:S467.
2. Haker S, Angenent S, Tannenbaum A, Kikinis R, Sapiro G, Halle M. Conformal surface parameterization for texture mapping. IEEE Transactions on Visualization and Computer Graphics. 2000 April–June;6(2):181–189.
3. Eck M, DeRose T, Duchamp T, Hoppe H, Lounsbery M, Stuetzel W. Multiresolution analysis of arbitrary meshes. Computer Graphics (Proceedings of SIGGRAPH 95); August 1995.
4. Kanai T, Suzuki H, Kimura F. Three-dimensional geometric metamorphosis based on harmonic maps. The Visual Computer. 1998;14(4):166–176.
5. Desbrun M, Meyer M, Alliez P. Intrinsic parameterizations of surface meshes. Proceedings of Eurographics. 2002
6. Vichnevetsky R. Computer Methods for Partial Differential Equations. Vol. 1. Prentice Hall; New Jersy: 1981.
7. Fischl B, Sereno MI, Tootell RBH, Dale AM. High-resolution inter-subject averaging and a coordinate system for the cortical surface. Human Brain Mapping. 1999;8:272–284. [PubMed]
8. Sisodiya S, Free S, Fish D, Shorvon S. MRI-based surface area estimates in the normal adult human brain: evidence for structural organization. J Anat. 1996;188:425–438. [PubMed]
9. Bakircioglu M, Grenander U, Khaneja N, Miller MI. Curve matching on brain surfaces using Frenet distances. Human Brain Mapping. 1998;6:329–333. [PubMed]
10. Thompson P, Toga A. A surface-based technique for warping 3-dimensional images of the brain. IEEE Transactions on Medical Imaging. 1996;15(4):1–16. [PubMed]
11. Shattuck D, Leahy R. Brainsuit: an automated cortical surface identification tool. Med Image Anal. 2002 Jun;6(2):129–42. [PubMed]
12. Gu X, Wang Y, Chan T, Thompson P, Yau S-T. Genus zero surface conformal mapping and its application to brain surface mapping. Information Processing in Medical Imaging. 2003 [PubMed]
13. Eells J, Sampson JH. Harmonic mappings of Riemannian manifolds. Ann J Math. 1964:109–160.
14. Woods R, Grafton S, Holmes C, Cherry S, Mazziotta J. Automated image registration: I. General methods and intrasubject, intramodality validation. Journal of Computer Assisted Tomography. 1998;22:139–152. [PubMed]
15. Thompson P, et al. Detecting Dynamic and Genetic Effects on Brain Structure using High-Dimensional Cortical Pattern Matching. Proceedings of ISBI. 2002 [PMC free article] [PubMed]