\documentclass{article}

\newcommand{\ba}{\left( \begin{array}}
\newcommand{\ea}{\end{array} \right)}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}
\newcommand{\bqn}{\begin{equation}}
\newcommand{\eqn}{\end{equation}}

\newcommand{\ot}{\displaystyle \circ}

\newcommand {\diff}[1] {\frac{\partial}{\partial #1}}
\newcommand{\jacob}[3]{\frac{\partial^2 #3}{\partial #1 \partial #2}}
\newcommand{\der}[2]{\frac{\partial #2}{\partial #1}}
\newcommand{\p}[1]{\partial_{#1}}

\renewcommand{\baselinestretch}{1.5}

\begin{document}
\newtheorem{lem}{Lemma}
\newtheorem{thm}{Theorem}
\newtheorem{defn}{Definition}

\title{\bf Shape Models of the Cortical Surface\\
\small Project for Prof. Collins' Advanced Medical Imaging Course}

\author{Moo K. Chung \\
\it{Department of Mathematics and Statistics}\\
\it{McGill University}\\
E-mail: chung@math.mcgill.ca \\}
\date{April 17, 2000}
\maketitle
\large
\begin{abstract}
We will examine the mathematical problem of how to represent the outer
cortex of the human brain. The cortical surface plays a important role
in several applications, including morphological analysis of the brain,
and spatial normalization and registration of brain images. Methods for
obtaining the Gaussian and the mean curvature will be presented in detail
which may be further used for extracting sulcal lines.  Then a statistical
analysis on the cortex can be done using the {\it Riemannian metric tensor}.
New method for
statistical analysis on local surface area change under nonlinear deformation
will be presented.
\end{abstract}

\section{Introduction}
The cerebral cortex, which makes up the largest part of the human
brain, has the topology of a 2D sheet, and a highly folded
geometry. Most of the features that distinguish these cortical
areas can only be measured relative to that local orientation of
the cortical surface. For this reason, surface-based methods is
desirable. However, most methods analyzing human brain imaging
data have relied upon 3D volume based approaches. For instance,
functional activation are typically reported in 3D volume of
human brain instead of on 2D cortical surface. But the amount of
buried cortex is estimated to be 60-70 percent, implying that
distances measured in 3D space between two points on the cortical
surface will substantially underestimate the true distance along
the cortical sheet, particularly in cases where the points lie on
different banks of a sulcus. The shortest distance in 2D or 3D
Euclidean space is a straight line while the shortest distance on
the cortical surface is a curved line called {\it geodesic
curve}. Therefore, it might be of certain interest to be able to
identify the geodesic curves between two different regions on the
cortex to see if functional connectivity is somehow related to the
distance measure on the surface. Such surface-based analysis can be
easily done if one can develop shape models on the cortex using
Riemannian metric tensor. \\

However, the major impediment to the routine use of cortical surface-based
approaches in brain imaging has been the difficulty of automating
the surface extraction, registration and analysis. In order to accomplish
the statistical analysis on the cortical surface, shape representation
of the cortical surface is an essential part. A cortical shape can be
represented in many ways. Triangular or polygonal mesh to form a wire frame
of the surface or continuous parameterization of the surface by certain
mathematical functions. In this paper, our emphasize is on the later
parametric representation of the curve and once such representation is
done how we approach the final task of the statistical analysis.
To facilitate this, we model the cortical
surface as a smooth two dimensional
manifold parameterized by $u$ and $v$ such that $$X(u,v) =
\{x_1(u,v),x_2(u,v),x_3(u,v):(u,v) \in D \}$$ represents the
parameterized surface on an open set $D \in \Re^2$. Before going into the
problem of parametric representation of the cortical surface,
consider the problem of parameterizing curves in $\Re^3$. Then as an extension
of the parameterization of curves, we can also parameterize surfaces.
The most widely used method of parameterizing an arbitrary curve in $\Re^3$
is called {\it B-spline interpolation}.

\section{B-Spline Interpolation of Curves and Surfaces}
A special B-Spline called {\it Bezier curves} was developed by Bezier
and De Casteljau(Gu\'{e}ziec,
1996) in the 1970's for CAD/CAM operations for car manufacturing industry.
It became the foundation
of the entire Adobe Postscript drawing model. We
can construct a parametric curve interpolating given points $V_0,\dots, V_m$
using Bezier curves. The simplest Bezier curve is a line segment
$V_{0,1}(\lambda)$ joining two points $V_0$ and $V_1$ and
parameterized by $\lambda \in (0,1)$
 such that
 $$V_{0,1}(\lambda) = (1-\lambda)V_0 + \lambda V_1$$
For three points $V_0,V_1$ and $V_2$, we
have a quadratic Bezier curve
$$V_{0,2}(\lambda) = (1-\lambda)^2V_0 + 2(1-\lambda)V_1 + \lambda^2 V_2.$$
In general, for $m+1$ points $V_1,\dots,V_m$, $m$-th order Bezier curve is
$$V_{0,m}(\lambda) = \sum_{j=0}^m
{m \choose j} (1-\lambda)^{m-j} \lambda^jV_j$$ Note that the basis functions
$\{ {m \choose j} (1-\lambda)^{m-j} \lambda^j,j=0,\dots,m \}$ are
{\it Bernstein polynomials} of degree $m$. Cubic Bezier curves are the
most often used Bezier curve because quadratic are necessarily
planner whereas cubics are not and cubics satisfy the minimum
curvature property(Gu\'{e}ziec, 1996) or minimum strain energy(McLeod, 1998)
which make them a more suitable
tool for data approximation. To avoid increasing the degree of
the Bezier curve, we need to piece together Bezier curves. If
continuity conditions are satisfied for each Bezier curve
segments, the result is a B-spline curve. In general, a B-spline
curve of degree $K-1$ with $m+1$ vertices $V_0,\dots V_m$ is
defined as $$S(u) = \sum_{j=0}^m B_{j,K}(u)V_j,$$ where the
B-spline functions $\{B_{j,k}(u),j=0,\dots, m \}$ are defined
recursively with $m+K+1$ knot values $u_j$ such that
\bq B_{j,1}(u) &=& 1, u_j \leq u < u_{j+1} \\
               & & 0, \mbox{ otherwise}
\eq
$$B_{j,k} = \frac{u-u_j}{u_{j+k-1} - u_j}B_{j,k-1}(u) +
\frac{u_{j+k}-u}{u_{j+k} - u_{j+1}}B_{j+1,k-1}(u).$$
For knot values $\{0,0,0,1,1,1\}$, we have the quadratic Bezier curve and
for knot values $\{0,0,0,0,1,1,1,1\}$, we have the cubic Bezier curve. \\

Then we can use these B-spline curves to generate surfaces.
Consider a rectangular mesh of vertices $\{V_{ij}\}$. A B-spline
patch $X$ of degree $K-1$ can be defined by the tensor product:
$$X(u,v) = \sum_{i,j}B_{i,K}(u)B_{j,K}(v)V_{ij}.$$
The advantage for using spline to represent the surface is that it is
easy to evaluate the curvature of a surface or other geometric
characteristics of the surface because polynomial functions can
be differentiated easily.

\section{Polynomial Regression of Surfaces}
Instead of using B-splines to form a parametric surface, there
is another method called {\it polynomial regression} (Ripple, 1981).
This is a smoothing technique to fit the given vertices $V_0,\dots,V_m$
by least squares method to a polynomial function of the form
$$f(x,y) = \sum_{i+j \leq p} a_{ij}x^iy^j$$

Then for $V_k=(x_k,y_k,z_k) \in \Re^3$, $P=\frac{(p+1)(p+2)}{2}$
unknown coefficients $a_{ij}$'s are chosen to minimize the residual
$$\sum_{k=0}^m [z_k - f(x_k,y_k)]^2.$$ The drawback of the polynomial
regression is that there is tendency to wave the outer most vertices
to find vertices in the center. So if we are warping the whole cortical
surface using a single polynomial function $f(x,y)$, the fit may break
down near the extreme points in $D$. However, from the purely statistical
point of view, we do not need to warp the cortex by a single
function $f$ over all domain $D$. What we really need is the parametric
representation in the neighborhood of each voxel. Later we will show that
the quadratic polynomial regression at each voxel with its
neighboring voxels is sufficient for analyzing the variability in cortex using
the Riemannian metric tensor approach. \\

Let $$f_1(X)=1,f_2(X)=x,f_3(X)=y,f_4(X) = x^2,$$
$$f_5(X)=xy,f_6(X)=y^2,\dots,f_P(X)=y^P.$$ and the coefficients
$a_{ij}$'s as $\beta_0,\beta_1,\beta_2,\dots, \beta_P$. Then the
problem becomes a standard {\it multiple regression analysis}:
$$z=X\beta + \epsilon,$$ where the matrix $X=(X_{ij})
=(f_j(x_i,y_i))$, the
parameter $\beta = (\beta_1,\dots,\beta_P)$ and $z=(z_1,\dots,z_k)$.
Then it could be given to any statistical package to estimate the
parameter $\beta$ automatically. Orthogonal polynomials are
often used for polynomial regression, especially with regularly spaced
data but modern multiple regression techniques rely on orthogonalizing
the functions $f_1,\dots,fp$ and make the orthogonal polynomial
technique unnecessary. Other families of surfaces has been suggested
but are rarely used. One of them is to use finite Fourier series: \bq
f(x,y) &=& \sum_{0 \leq i,j \leq r}
[a_{ij}\sin(i\omega_1x)\sin(j\omega_2y)
+ b_{ij}\sin(i\omega_1x)\cos(j\omega_2y) \\
&& + c_{ij}\cos(i\omega_ix)\sin(j\omega_2y)
+ d_{ij}\cos(i\omega_1x)\cos(j\omega_2y)]
\eq
Such a surface requires a large number of coefficients,
i.e. $4r^2 +4r+1$, so even with $r=3$ there are 49 coefficients
plus the fundamental frequencies $\omega_1$ and $\omega_2$ to be
estimated. Hence, surface warping based on Fourier series on the large
data set is computationally intensive although the result might be better than
the polynomial regression. With small values of $r$, the surface shows ripples at
the fundamental frequencies in the $x$ and $y$ directions and so is
really suitable for periodic surface. In certain sense, the cortical
surface is somewhat periodic in a sense that there is a continuous ups
and downs of the cortical folds. So the above Fourier series expansion can be
used to represent the cortex.
Instead of using the flat domain $D \in \Re^2$, we can also
find a smooth mapping from the cortical surface to a unit sphere
using the spherical harmonics,
which allows points on the sphere to be represented by two angles
$\theta$ and $\psi$ (Bakircioglu {it et al.}, Thompson {\it et al.}, 1996)

\section{Quadratic Surfaces}
Define the tangent plane $T_p$ at a point $p=(x,y,z)$ on the
cortical surface having the normal vector $n_p$. Let $b_1$ and
$b_2$ be a orthonormal basis of the tangent plane. The quadratic
surface approximating the parametric surface passing thru the
point $p$ can be written as $$ X(p) = p + ub_1 + vb_2 +
[u,v]C_p[u,v]^t n_p,$$ where $C_p$ is a $2 \times 2$ symmetric
matrix. Note that $X(p)$ is a quadratic polynomial in $u$ and $v$.
For any parametric surface $X = X(u,v)$, most differential geometric
properties can be represented using up to the second partial derivatives
$X_{uu},X_{uv}$ and $X_{vv}$. Therefore, the quadratic polynomial regression
on each voxel and its neighboring voxels is essentially all we need to
analyze the geometry of the cortical surface. \\

If a minimum mean squared error
estimate in the neighborhood of the point $p$ is used, then $C_p$
becomes a curvature matrix and the {\it Gaussian curvature} $K_G$
and the {\it mean curvature} $K_H$ at the point $p$ is given by
$$K_G = det(C_p), K_H = tr(C_p)/2.$$
 The {\it principal curvatures} are defined as the maximum and
 minimum cross-sectional curvature of the surface and they are
 given by the eigenvalues of the curvature matrix $C_p$. Hence,
 the Gaussian curvature is the product of the principal curvatures
 while the mean curvature is the mean of the principal curvatures.
 The estimation of the Gaussian and the mean curvature can be numerically
 computed on a triangulated mesh representing the cortical surface by
 Hamann's scheme(Joshi {\it et al.}).

\section{Differential Geometry on the Cortical Surface}
The surface of the outer cortex may be assumed
to be a two dimensional manifold of class $C^2$, i.e.
$$X(u,v) = \{ x_1(u,v),x_2(u,v),x_3(u,v): (u,v) \in D \subset \Re^2 \}$$
such that all partial derivatives of order two or less are
continuous in $D$. The pair of the mapping $X$ and its domain $D$
is called a coordinate neighborhood or a {\it local chart}.
A local chart can be approximated at each voxel in the surface by
locally approximating the diffeomorphism by up to the quadratic
terms in its Taylor series using the above quadratic polynomial regression.
If $D$ is a unit square in $\Re^2$ and the surface $S$ is smoothly deformable
from a sphere then at least two different local charts are required.
The partial derivative vectors $X_u$ and $X_v$
of the surface $X(u,v)$ with respect to $u$ and $v$ and form the
{\it covariant} base vectors on the tangent plane. If $$X_u
\times X_v \neq 0$$ in a region on the surface, a single-valued
inverse of the transformation exists in that region. For
convenience, let $(u^1,u^2)=(u,v)$. Then the {\it Riemannian
metric tensor} $g_{ij}$ on the surface is given by $g_{ij} =
<X_i, X_j>,$ where the indices $i,j$ are either $u$ or $v$. Note
that $G=(g_{ij})$ forms $2 \times 2$ symmetric matrix. Let $g = \det(G)=
g_{11}g_{22} - g_{12}^2$. The unit normal vector $n$ to the
surface is $$n = \frac{X_u \times X_v}{\sqrt{g}}.$$ Then the
vector $(X_u, X_v, n)$ forms a right-handed triad in $\Re^3$. The
curvature tensor $K_{ij}$ is given by $$K_{ij} = \der{u^j}{X_i}
\cdot n.$$ Note that $K=(K_{ij})$ forms the curvature matrix and
the Gaussian and the mean curvature are given by the determinant
and the trace of $K$ as in the case of $C_p$. The metric tensor
and the quadratic surface approximation approaches are
equivalent(Joshi {\it et al.}). The metric tensor method is very general
and can be used for computing not only the curvatures but also
other geometric properties such as the surface areas or geodesic curves
on the cortex.

\section{Estimating the Riemannian Metric Tensor}
If we have a parametric representation $X=X(u,v)$ of the cortical surface,
calculating the metric tensor is trivial. However, the parametric
representation of the surface itself is the hardest part of the problem. One
approach to avoid the difficulty is not to parameterize the cortical
surface by continuous functions directly but to warp the surface with a fine
triangular or polygonal mesh first and then use this triangular mesh to find
the parametric representation(Joshi {\it et al.}). The standard method for triangulating the
surface is the {\it marching cubes} algorithm. Alternating methods such as
Coordinate Curves(Timsari), Deformable Surface Algorithm(Davatzikos, 1996)
and Double Hull Algorithm(MacDonald, 1997) are available. \\

However, if estimating the metric tensor is the main problem, there
is a new method without using any parameterization.
Let $S=\{(x,y,z,w): w \in \{0,1\} \}$ be the extracted
cortex from a 3D brain image with the binary mask $w \in \{0,1\}$. So at each
voxel $p=(x,y,z)$, if $p$ is the part of the cortex then the voxel has
a value $1$ and $0$ otherwise. Let $D_z$ be the projection of the $S$
onto the plane $z=0$. Then the trivial mapping
$$X(x,y) = \{ x(x,y)=x, y(x,y)=y, z(x,y)=z: (x,y) \in D_z \}$
with the domain $D_z$ forms a local chart if $X$ is a $C^2$ function.
Similarly, we can form
local chats $D_x$ and $D_y$ by projecting $S$ onto the plane $x=0$ and $y=0$
 respectively.
Assuming the cortical surface $S$ to be deformable to a unit sphere
by a smooth mapping, at least two local charts are required. Hence,
three local charts $D_x,D_y,D_z$ are more than adequate to cover the whole
cortex. Using the local chart $D_z$, the covariant base vectors
are given by $X_x = (1,0,z_x)'$ and $X_y =(0,1,z_y)'$.
So the metric tensor is
$$g_{11}=1+z_x^2, g_{12} = z_xz_y, g_{22} = 1 + z_y^2$$
Then we can show that the Gaussian curvature is given by
$$ K = \frac{z_{xx}z_{yy} - z_{xy}^2}{(1+z_x^2+z_y^2)^2}.$$
The local chart $D_z$ breaks down when $z_x = \infty$ or $z_y=\infty$. If
one of the partial derivatives $z_x$ or $z_y$ goes to $\infty$
then simply choose different local chart $D_y$ or $D_z$. This happens
when the tangent plane on the cortex is vertical. For $C^2$ manifold $S$,
I claim that the three local charts $D_x,D_y,D_z$ can not break down at the same
time and at least one of the local charts can be used.

\section{Application 1: Geodesic Curves}
The most major application of the metric tensor approach so far in the literature
has been an extraction of sulcal lines by maximizing the Gaussian
curvature(Bakircioglu) and a cortical surface registration(Thompson, 1996).
As another possible application of the metric tensor approach we can analyze
the arc length of an arbitrary curve $C$ lying on the
cortical surface such as sulcal lines. Parameterize the curve $C$
by the arc length $s$ measured along the curve from some fixed
point on $C$, and suppose that the surface coordinates of $C$ can
be expressed parametrically in terms of $s$: $$x^1=x^1(s),
x^2=x^2(s) \mbox{ along } C.$$ Then the curve of shortest
distance joining two points on a surface called a {\it geodesic}
is given by the differential equation: $$\frac{d^2
x^k}{d s^2} + \Gamma_{ij}^k \frac{d x^i}{ds}\frac{d x^j}{ds} =
0,$$ where the Christoffel symbol with $G^{-1} = (g^{ij})$ is given by,
$$\Gamma_{ij}^k = \frac{1}{2}g^{ks}(\der{x^j}{g_{si}}-\der{x^s}{g_{ij}} +
\der{x^i}{g_{js}}).$$
The equation can be solved numerically so that it is possible to
find a path connecting one region of the cortical surface to another
by the shortest path and it would be interesting to see if there is any
correlation with a functional connectivity on the cortex.

\section{Application 2: Local Surface Area}
The study of statistical variability in local area change in the
2D cross sectional image of the corpus callosum has been done by
Davatzikos {\it et al.} in investigating sex differences in the
corpus callosum. The deformation is from 2D Euclidean plane to the
target 2D Euclidean plane. So in the case of the cortical surface,
Davatzikos' method does not work and we need new mathematical treatment
based on the metric tensor of the curved manifold. As brain grows over time,
the outer cortical surface itself expand as well as the interior whole brain
volume. It is highly likely that the cortical surface expansion or
shrinking is not uniform. Also in comparing male and female cortical
surface, certain region of the cortical surface of male brains might
be smaller or larger than that of female brain. One measure of testing such
variability is the {\it local surface area}.\\

The element of the surface area $dS$ is given
by(Kreyszig, 1964) $$dS = \sqrt{g}dudv.$$ Integrating over $D$,
we get the total cortical surface $S$, i.e. $$S = \int_D
\sqrt{g}dudv.$$ Note that $\sqrt{g}$ measures the area of the
unit square after the transformation from $D$ to $S$ by $$(u,v)
\to X(u,v).$$ If the surface $X(u,v)$ deforms to $X(u,v) +
U(X(u,v)),$ where $U$ is the displacement vector field, then the
local surface area $\sqrt{g}$ also changes accordingly. To
measure the change in the local surface area, we need to know how
the metric tensor $g_{ij}$ changes by the deformation. Let
$Y(u,v)$ be the deformed surface by the displacement vector field
$U$, i.e.
$$Y(u,v) = X(u,v) + U(X(u,v)).$$
Then the partial derivative vector is
$$Y_i = X_i + \nabla U \cdot X_i,(i=u,v)$$
where $\nabla U = (\der{x_j}{U_i})$. Hence, the metric tensor
$h_{ij}$ on the deformed surface $Y(u,v)$ is
\bq h_{ij} &=& <Y_i,Y_j> \\
           &=& <X_i,X_j> + 2X_i'\nabla U X_j + X_i'(\nabla U)' \nabla U X_j \\
&=& g_{ij} + 2X_i' \nabla U X_j + O(\der{x_k}{U_l}\der{x_m}{U_n}).\eq
When there is no deformation, $h_{ij} = g_{ij}$. Let $H=(h_{ij})$ and $h= det(H) = h_{11}h_{22} - h_{12}^2$. For small displacement, the metric tensors in matrix formulation is given by
$$H \approx G + 2(\nabla X)' \nabla U \nabla X,$$
where $\nabla X = (X_1, X_2)$.
Then, the change of local surface area element per unit surface area
denoted by $\lambda(u,v)$ is

\bq \lambda(u,v) = \frac{\sqrt{h}-\sqrt{g}}{\sqrt{g}} &=& \sqrt{\frac{h}{g}}-1 \\
&=& [\det(HG^{-1})]^{1/2}-1\\
&=& [\det( I + 2(\nabla X)'\nabla U \nabla X G^{-1})]^{1/2} -1\\
&\approx& 2tr[(\nabla X)' \nabla U \nabla X G^{-1}]
\eq
Assuming $U(u,v)$ to be 2D Gaussian random vector field on $D$,
$\lambda(u,v)$ is a Gaussian field and the detection on the local
surface area change can be performed the standard two sample $t$ test.

\section{Application 3: Cortical Surface Matching}
Once a parametric representation of the outer cortical surface has been
obtained, geometric features of the cortical surface, such as various
curvatures, can be determined and used for the identification and
matching of sulci, gyri, fissures, and other distinct cortical
characteristics. Then the surface deformation derives 3D volumetric
warping of the rest of the brain(Davatzikos, 1996, Thompson and Toga, 1996).
3D active surfaces(Cohen and Cohen, 1992) are used to automatically extract
parametric representations of each subject's cortex, on which
corresponding networks of anatomical curves are identified.
The transformation relating these networks is expressed as a vector
field in the parameter space of the cortex. This vector field in parameter
 space indirectly specifies a correspondence field in 3D, which drives
 one cortical surface into the shape of another. This mapping matches up
 every element of the specified system of landmark curves with its
 counterpart in the target brain.


\section{Conclusions and Future Direction}
By modeling the cortical surface as two dimensional smooth manifold,
the Riemannian metric tensor can be computed on the surface which can be further
used for classification, registration and statistical analysis. My emphasis on
this paper has been to present different methods for a parameterization the cortex
using a given triangular mesh, how to compute the metric tensor and
its applications. While investigating the field, I have come up with
two new ideas which have been presented in section 6 and section 8.
In the future, I hope to actually compute
the metric tensor numerically using the method in section 6 and
apply the result for testing inter- and intra-subject variability of the cortex. \\



{\bf Note on References}: {\it Some of the papers are from the web
page of the authors, so they may not have proper year and journal
where they are published.}
\begin{thebibliography} \\
\bibitem{Bakircioglu} MM Bakircioglu, S Joshi, MI Miller, {\sl Landmark Matching on Brain Surfaces via Large Deformation Diffeomorphisms on the Sphere}
\bibitem{Dale} AM Dale, B Fischl, {\sl Cortical Surface-Based Analysis I: Segmentation and Surface Reconstruction.}
\bibitem{Fischl} B Fischl, MI Sereno, {\sl Cortical Surface-Based Analysis II: Inflation, Flattening, and a Surface-Based Coordinate System.}
\bibitem{Davatzikos} C Davatzikos RN Bryan, {\sl Using a Deformable Surface Model to Obtain a Shape Representation of the Cortex}, IEEE Trans. Medical Imaging, Vol. 15, No. 6, 1996
\bibitem{Davatzikos2} C Davatzikos, {\sl Nonlinear Registration of Brain Images Using Deformable Models}, IEEE, 1996
\bibitem{Davatzikos3} C Davatzikos, M Vaillant, SM Resnick, JL Prince, S Letovsk, RN Bryan,{\sl A Computerized Approach for Morphological Analysis o the Corpus Callosum}
\bibitem{Gueziec} A Gu\'{e}ziec, {\sl Spline Curves and Surfaces for Data Modeling}, Advances in Morphometrics, Plenum Press, New Work, 1996
\bibitem{Joshi} SC Joshi, J Wang, MI Miller, {\sl On the differential geometry of the cortical surface}
\bibitem{Khaneja} N Khaneja, MI Miller, U. Grenander {\sl Dynamic Programming Generation of Curves on Brain Surfaces}
\bibitem{Ripley} BD Ripley, {\sl Spatial Statistics}, John Wiley $\&$ Sons, Inc, 1981
\bibitem{Kreyszig} E Kreyszig, {\sl Differential Geometry}. University of Toronto Press, Toronto, 1964
\bibitem{MacDonald} JD MacDonald, {\sl A Method for Identifying Geometrically Simple Surfaces from Three Dimensional Images}, Ph.D. Thesis, McGill University, 1997
\bibitem{McLeod} RJY McLeod, ML Baart, {\sl Geometry and Interpolation of Curves and Surfaces}, Cambridge University Press, 1998
\bibitem{Thompson} P Thompson, AW Toga, {\sl A Surface-Based Technique for Warping Three-Dimensional Images of the Brain}, 1996

\end{thebibliography}
\end{document}