Xiuwen Liu

Washington Mio

Yonggang Shi

Ivo Dinov

Geodesic Shape Spaces of Surfaces of Genus Zero ⋆

We construct shape spaces of elastic spherical surfaces immersed in Euclidean space R k . The spaces are equipped with geodesic metrics that depend on the tension and rigidity of the surfaces. We develop algorithms to calculate geodesics and geodesic distances, as well as tools to quantify local shape similarities and contrasts, thus obtaining a local-global formulation. We give examples of geodesic interpolations and illustrations of the use of the model in brain mapping.

Introduction

The development of computational models of shapes of surfaces in Euclidean space is motivated by core problems in computational anatomy such as mapping the human brain, characterizing normal variation in anatomy and identifying pathological changes. In recent years, significant progress has been made in the study of Riemannian shape spaces of curves (see e.g. [START_REF] Younes | Computable elastic distance between shapes[END_REF][START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF][START_REF] Klassen | Analysis of planar shapes using geodesic paths on shape manifolds[END_REF][START_REF] Mio | On shape of plane elastic curves[END_REF][START_REF] Michor | A metric on shape space with explicit geodesics[END_REF]). However, a computationally feasible approach to higher-dimensional shapes still poses many challenges. In this paper, we construct a geodesic shape space of elastic surfaces of genus zero in R k , with shapes realized as piecewise linear immersions of a triangle mesh K homeomorphic to the 2-sphere. We represent a parametric shape by its discrete exterior derivative because first-order representations provide a good balance between computational tractability and geometrical accuracy. Representations beyond order zero also are more effective in capturing deformations such as elastic bends and folds. The space of immersions is equipped with a family of Riemannian metrics that reflect the elasticity properties of the surfaces. The general model is anisotropic and inhomogeneous, as it allows the tension and rigidity of a surface to vary throughout its extension. The representation is invariant under translations and a shape space of spherical surfaces is obtained by normalizing scale and accounting for the action of orthogonal transformations.

The paper is devoted to the development of the elastic shape models and algorithms to calculate geodesic paths and geodesic distances in shape space. We use these basic tools to estimate mean shapes and to construct an anatomical atlas of the right hippocampus. Our shape model is used in conjunction with existing surface parametrization and registration techniques [START_REF] Praun | Spherical parametrization and remeshing[END_REF][START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF][START_REF] Shi | Direct mapping of hippocampal surfaces with intrinsic shape context[END_REF], but research in these areas is still ongoing and is complementary to the study of shape spaces. Although shape distance is a global measurement of shape dissimilarity, the formulation used also enables us to characterize and quantify local shape differences. This is important in applications in order to identify the regions where the most significant morphological contrasts occur. Shapes of surfaces can be studied from several other viewpoints. Alternative approaches include: (i) models based on landmark representations [START_REF] Kendall | Shape manifolds, Procrustean metrics and complex projective spaces[END_REF]; (ii) models that rely on diffeomorphisms acting on a volume containing a surface to describe shape deformations [START_REF] Grenander | Computational anatomy: An emerging discipline[END_REF]; (iii) models based on medial-axis representations [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF].

The paper is organized as follows. In Section 2, we discuss the shape representation and examine invariance under shape-preserving transformations. In Section 3, we define the elastic shape metrics. The algorithm to calculate geodesics via energy minimization is presented in Sections 4 and 5, where illustrations are also provided. Energy density functions and localization techniques are discussed in Section 6. We conclude with a summary and some discussion in Section 7.

Shape Representation

Let K be a finite simplicial complex whose underlying polyhedron |K| is homeomorphic to the unit 2-sphere S 2 . K will be fixed throughout the discussion. A parametric spherical shape will be realized as a mapping α : |K| → R k , which is linear on each simplex of K. Thus, if V = {z 1 , . . . , z p } is the vertex set of K, α is completely determined by its restriction α : V → R k to V . Next, we define the discrete exterior derivative dα, which is a measure of the variation of α along the edges of K. Fix an orientation for each edge and each triangle of K and let E = {e 1 , . . . , e m } and F = {T 1 , . . . , T n } be the oriented edge and face sets of K. For an oriented edge e, let e -and e + denote its initial and terminal vertices, respectively. Then, dα : E → R k is given by dα(e i ) = α(e + i ) -α(e - i ). Note that it suffices to define dα over the oriented edges in E because if we reverse the orientation of an edge, the variation of α gets multiplied by -1. We only consider mappings such that dα(e i ) = 0, for every e i ∈ E, which we refer to as immersions. This just means that no edge gets crushed to a point under α. For each i, 1 i m, we write the modular component of dα(e i ) in logarithmic scale as r i = log dα(e i ) and the directional component as v i = dα(e i )/ dα(e i ) . Thus, the (oriented) ith edge of α is represented by the vector e ri v i . Writing each v i as a row vector v i = [v i1 . . . v ik ], the parametric surface α will be represented by the pair (r, v) ∈ R m × R m×k , where r = r 1 . . . r m T and v is the m × k matrix whose ith row is v i . Since the Riemannian metrics on R m × R m×k to be used in the development of our shape models differ from the standard Euclidean metric, we will use the notation E for this space to emphasize this fact. Thus, E is the space of dimension m(k + 1) formed by all pairs (r, v), with the topology induced by the Euclidean metric. Thus far, the only constraint on (r, v) is that each row of v must be a unit vector. We denote by M the subspace of E formed by all pairs (r, v) with this property.

Pre-Shape Space

The representation of α via (r, v) is clearly invariant under translations. To fix scale, we set the total edge length to be unitary; that is, we normalize the pair (r, v) to satisfy

G 1 (r, v) = m i=1 e ri = 1. (1)
Remark. Alternatively, one may normalize scale by fixing the total area, or simply drop this condition if a scale-sensitive model is desired. What pairs (r, v) represent the exterior derivative of a mapping α : |K| → R k ? The integrability conditions will reveal the further constraints to be imposed on (r, v). If p is an oriented path in K formed by a sequence of oriented edges, let the integer n i be the net number of times that the oriented edge e i ∈ E is traversed by p, where a negative sign indicates reversal of orientation. Given an immersion α, its variation along p may be expressed as m i=1 n i dα(e i ) = m i=1 n i e ri v i . The variation of α along any (oriented) cycle c in K clearly must vanish. Conversely, given (r, v), define the integral of (r, v) along the path p to be I(p) = i n i e ri v i . If I(c) vanishes along every cycle c, then (r, v) represents the exterior derivative of some α, uniquely determined up to translations. This can be seen as follows: fix a vertex v 0 of K and a point x 0 ∈ R k and define α(v 0 ) = x 0 . For any other vertex v, choose a path p from v 0 to v and let α(v) = x 0 + I(p). The vanishing condition over cycles ensures that α(v) is independent of the path chosen.

Verifying that I(c) = 0, for every cycle c, is not computationally feasible. However, it suffices to check this vanishing condition on cycles that are boundaries of oriented triangles (this uses the fact that the 2-sphere is simply connected). This is because the integral along any cycle in K can be accounted for by combining the integrals along the boundaries of appropriately chosen oriented triangles. For each oriented triangle T ℓ ∈ F , let ε ℓ1 e ℓ1 , ε ℓ2 e ℓ2 , ε ℓ3 e ℓ3 be the oriented edges of T ℓ , where ε ℓi = ±1, e ℓi ∈ E, and the indexes satisfy ℓ 1 < ℓ 2 < ℓ 3 . Then, the argument above shows that a pair (r, v) represents the exterior derivative of an immersion α if and only if it satisfies the following kn integrability conditions:

G ℓ,j (r, v) = 3 i=1 ε ℓi e r ℓ i v ℓi,j = 0 . (2)
where 1 ℓ n and 1 j k. It is not difficult to see that one may drop condition (2) over one of the triangles, say T n , since the variation along the boundary of any of the triangles can be expressed as a combination of the variations along the boundaries of the remaining ones. Thus, we obtain k(n -1) independent integrability conditions. Pairs (r, v) ∈ M satisfying (1) and (2) will be called pre-shapes and the space of all pre-shapes will be denoted P . A pre-shape gives a representation of immersions that is invariant to scale and translations. Thus, we have introduced a nested sequence P ⊂ M ⊂ E of spaces, which will be useful at various stages of our constructions.

Remark. A more formal explanation of the integrability of (r, v) can be given in terms of simplicial cohomology. The vanishing of the integral of (r, v) along the boundaries of oriented triangles means that the simplicial 1-cochain (with coefficients in R k ) determined by e i → e ri v i , 1 i m, is a 1-cocycle. This cocyle represents the trivial cohomology class because |K| is simply connected. Thus, the cocycle is integrable. This interpretation also allows us to identify the further integrability conditions needed for a complex K of different topology. The details will be discussed in future work.

Shape Space

Let O(k) be the group of k × k orthogonal matrices. If U ∈ O(k) and α is an immersion, the result of applying the rigid transformation

U to α is the composition U • α : V → R k . The induced action on (r, v) is (r, v) → (r, vU T ).
The action on the modular component is trivial because edge lengths do not change under U . Pre-shapes clearly get transformed into pre-shapes and the orbit of (r

, v) is O(r, v) = {(r, vU T ) : U ∈ O(k)}.
Since pre-shapes that differ by a rigid transformation have the same shape, we define the shape space S as the orbit space of P under the action of O(k). In other words, as the quotient S = P/O(k).

Riemannian Metrics

We introduce Riemannian structures on E that will lead to a family of elastic preshape and shape metrics. The general model is anisotropic and inhomogeneous, but the family contains a special 1-parameter homogeneous class. The idea is that each edge of a shape may offer different resistance to stretching/compression and bending, which will be quantified by its tension and rigidity coefficients. The geodesic distance will be related to the minimum energy required to deform a shape into another under these conditions.

Let a, b : E → R + be positive functions defined on the edge set, where a i = a(e i ) and b i = b(e i ) represent the tension and rigidity coefficients of the ith edge. Consider the Riemannian structure on E given at (r, v) by

(h, w), ( h, w) (r,v) = m i=1 a i h i hi e ri + m i=1 b i (w i • wi )e ri . (3)
If we think of (h, w) as an infinitesimal deformation of (r, v), then (h, w) 2 (r,v) = (h, w), (h, w) (r,v) is the energy cost of the deformation. Henceforth, the submanifolds P and M will have the Riemannian structure induced from (3). If (r, v) and (r, v) are pre-shapes, we let δ((r, v), (r, v)) be the geodesic distance between them on the pre-shape manifold P . If s and s are the shapes associated with the orbits of (r, v) and (r, v), the shape distance is defined as

d(s, s) = inf U,V ∈O(k) δ((r, vV T ), (r, vU T )) = inf U ∈O(k) δ((r, v), (r, vU T )) . (4)
The last equality follows from the fact that the Riemannian metric is compatible with the action of O(k); that is, O(k) acts by isometries. Most of the remaining work will be devoted to the calculation of geodesics in P and the geodesic distance δ. The minimization over O(k) is a relatively simpler calculation, as indicated below.

Path Spaces and the Energy Functional

Our next goal is to construct a pre-shape geodesic between (r, v), (r, v) ∈ P . The idea is to begin with a path in P connecting the pre-shapes and gradually deform it to a geodesic following the negative gradient flow of the path energy. Implicit in this statement is that a path space equipped with a Riemannian structure has been constructed.

Path Spaces

Let I = [0, 1]. If p : I → M is a path in M , we write its modular and directional components as p(t) = (r t , v t ). Let Y be the manifold of all absolutely continuous paths in M (with square integrable derivative). A tangent vector to Y at p represents an infinitesimal deformation of the path and is given by a vector field (h t , w t ) along p with the property that w t i • v t i = 0, for every 1 i m and ∀t ∈ I. This follows from the fact that v t i 2 = 1 so that the deformation of each v t i has to be tangential to the unit sphere. Consider the Riemannian metric on the path space Y given by (h, w), ( h, w) p = (h 0 , w 0 ), ( h0 , w0 ) p(0)

+ 1 0 D t h t , w t , D t ht , wt p(t) dt , (5)
where D t denotes covariant differentiation in M along p. This type of inner product in function space was, to our knowledge, introduced by Palais [START_REF] Palais | Morse theory on Hilbert manifolds[END_REF]. Motivated by the inclusion P ⊂ M , we consider the following submanifolds of Y :

(i) The space Z M ⊂ Y of paths in M satisfying the boundary conditions p(0) = (r, v) and p(1) = (r, v). (ii) The space Z P ⊂ Z M of paths in the pre-shape space P satisfying the boundary conditions described in (i). A path in Z P has the property that each (r t , v t ) also satisfies (1) and (2).

We thus have three nested path spaces:

Z P ⊂ Z M ⊂ Y .

The Energy Functional

On the path space Y , define the energy functional E : Y → R by

E(p) = 1 2 1 0 (∂ t r t , ∂ t v t ), (∂ t r t , ∂ t v t ) p(t) dt . (6)
A pre-shape geodesic between (r, v), (r, v) ∈ P is a path p ∈ Z P , which is a critical point of the energy E restricted to Z P . We are particularly interested in minimal energy paths since they represent minimal-length geodesics. We exploit the inclusions Z P ⊂ Z M ⊂ Y in our approach to the minimization problem, the point being that it is easier to calculate the gradient of E as a functional on Z M . The use of ( 5) leads to a particularly elegant and computationally robust expression for the gradient.

Geodesics

In this section, we develop an algorithm to calculate a geodesic in P between two pre-shapes (r, v) and (r, v) using a gradient search to minimize the energy.

Initialization

To initialize the process, let α, ᾱ : V → R k be immersions of K associated with (r, v) and (r, v), respectively. We linearly interpolate α and ᾱ to obtain a 1parameter family of mappings α t : K → R k . If dα t (e i ) = 0, we gently deform α t to make it non-singular and then scale each α t to turn dα t into a pre-shape. The path p(t) = dα t is used to initialize the search.

Covariant Integration and Parallel Transport

To calculate the gradient of E at p, we first discuss covariant integration in M of a vector field (f t , x t ) along a path p(t) = (r t , v t ), where f t , where x t denote the modular and directional components of the field. The tangentiality to M of the field (f t , x t ) just means that x t i • v t i = 0, ∀t ∈ I and every i, 1 i m. As shown in Appendix B, a vector field (F t , X t ) along p is tangential to M and is a covariant integral of (f t , x t ) if and only if it satisfies the system of differential equations

∂ t F t i = f t i -1 2 (∂ t r t i )F t i + 1 2 bi ai (X t i • ∂ t v t i ) ∂ t X t i = x t i -1 2 (X t i ∂ t r t i + F t i ∂ t v t i ) -(X t i • ∂ t v t i )v t i . (7)
Numerical integration of (7) with initial conditions (F 0 , X 0 ) will be used in the calculation of geodesics. In the special case where the field (f t , x t ) is identically zero, the integral field (F t , X t ) is the parallel transport of (F 0 , X 0 ) along p.

The Gradient of the Energy

Given a path p ∈ Z P , we first calculate the gradient of E at p as a functional on the path space Y . For this purpose, we consider a variation (r t (µ), v t (µ)) of p(t) = (r t , v t ) in Y along a direction (h t , w t ), where µ ∈ (-ǫ, ǫ) is the variation parameter. This means that (r t (0), v t (0)) = (r t , v t ), the path (r t (µ), v t (µ)) ∈ Y for each fixed µ,

h t = ∂ ∂µ r t (µ) µ=0 and w t = ∂ ∂µ v t (µ) µ=0 . (8)
Differentiating ( 6) at µ = 0, we obtain

dE p (h t , w t ) = 1 0 D t (h t , w t ), (∂ t r t , ∂ t v t ) p(t) dt ,
where D t denotes covariant derivative. If we set (f t , x t ) = (∂ t r t , ∂ t v t ) in ( 7) and integrate the system with initial condition (F 0 , X 0 ) = (0, 0), we get a vector field (F t , X t ) along the path p. Then, we may rewrite the variation dE p (h t , w t ) as dE p (h t , w t ) = (h t , w t ), (F t , X t ) p .

Thus, the gradient of E at p as a functional on Y is ∇ Y E(p) = (F t , X t ). To obtain ∇ Z M E(p), we project ∇ Y E(p) orthogonally onto the tangent space of Z M at p with respect to the inner product [START_REF] Michor | A metric on shape space with explicit geodesics[END_REF]. Since the space Z M is obtained from Y by imposing the boundary conditions (r 0 , v 0 ) = (r, v) and (r 1 , v 1 ) = (r, v) on paths, a tangent vector to Y at p is tangent to Z M if and only if it vanishes at t = 0 and t = 1. A simple (covariant) integration by parts argument shows that the orthogonal complement of the tangent space of Z M in the tangent space of Y at p is formed by the covariantly linear fields (that is, smooth fields with trivial second covariant derivative) in M along p. By construction, the field (F t , X t ) vanishes at t = 0. Thus, to orthogonally project (F t , X t ) onto the tangent space of Z M , we simply need to subtract from (F t , X t ) the covariantly linear field that vanishes at t = 0 and coincides with (F t , X t ) at t = 1. Again, we resort to covariant integration, this time applied to the reverse of the path p. First, construct a covariant field by integrating the everywhere zero field with initial condition (F 1 , X 1 ) along the reverse of the path p. Reversing the path again, after integration, we obtain a parallel field (G t , Y t ) along p whose value at t = 1 is (F 1 , X 1 ). The field (tG t , tY t ) is covariantly linear with the desired properties. Therefore, the gradient is given by

∇ Z M E(p) = ∇ Y E(p) -(tG t , tY t ) . ( 10
)
Our goal is to minimize E on the path space Z P . One possible approach is to calculate ∇ Z P E(p), which would allow us to implement a gradient descent directly in Z P . However, in this case, the calculation is quadratic in the number m = |E| of edges. This is undesirable as spherical meshes used in neuroimaging often have a large number of edges. Thus, we employ an alternative numerical strategy that scales linearly with m and is based on the replacement of gradient descent in Z P by its counterpart in Z M followed by a closest-shape projection onto Z P . The projection algorithm is discussed in Appendix A.

Pre-shape Geodesics

We now present an algorithm to estimate a pre-shape geodesic from (r, v) to (r, v). Let ǫ, δ > 0 be small real numbers:

(i) Initialize the search with a path p(t) = (r t , v t ) in Z P , e.g., as described in Section 5. (ii) Let (f t , x t ) = (∂ t r t , ∂ t v t ). Using [START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF], integrate this field covariantly along the path p with zero initial condition. The integral field (F t , X t ) gives the gradient ∇ Y E(p). (iii) Using [START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF], calculate the parallel transport of (F 1 , X 1 ) along the reverse of the path p. Reverse the path and the parallel field again to obtain a parallel field (G t , Y t ) along p. By [START_REF] Grenander | Computational anatomy: An emerging discipline[END_REF], the

Z M -gradient of E is ∇ Z M E(p) = ∇ Y E(p) - (tG t , tY t ).
(iv) Write the modular and the directional components of the gradient

∇ Z M E(p) as (h t , w t ), 0 t 1. Update p(t) = (r t , v t ) as a path in Z M according to (a) rt = r t -ǫh t ; (b) ṽt i = v t i , if w t i = 0; (c) ṽt i = cos(ǫ w t i )v t i - sin(ǫ w t i ) w t i w t
i , otherwise. Note that the update of v t i takes place along great circles of the unit sphere in R k ensuring that each ṽt i is a unit vector. (v) Project each (r t , ṽt ) onto the pre-shape space P (see Appendix A) to obtain a path (r t , v t ) new in P . (vi) Iterate the process until (r t , v t ) new -(r t , v t ) p < δ, where • p denotes the Palais norm defined in (5).

Shape Geodesic and Shape Distance

The geodesic shape metric defined in (4) involves a minimization of the preshape distance over the orthogonal group O(k). We briefly indicate how the minimization problem can be treated. We obtain an initial estimate of Û = argmin

U ∈O(k) δ((r, v), (r, vU T )) , (11)
for example, by minimizing (r, v) -(r, vU T ) 2 (r,v) . This is the same as minimizing

m i=1 b i v i -vi U T 2 e ri , which is equivalent to maximizing the expression m i=1 (b i e ri/2 v i ) • (v i U T
). This problem is analogous to the one that arises in Procrustes alignment of shapes [START_REF] Kendall | Shape manifolds, Procrustean metrics and complex projective spaces[END_REF] and admits a closed form solution. Once this initial estimate is obtained, a gradient search can be used to locally refine the estimation. Using the algorithm to calculate pre-shape distances, we compute the O(k)-gradient numerically via finite differences using the Lie group structure of O(k).

Figures 1(a) and (b) show geodesics between cortical surfaces with 122,880 edges and hippocampal surfaces with 30,720 edges. The geodesic interpolations were calculated with homogeneous elasticity coefficients a i = 0.3, b i = 0.7 for the cortices and a i = 0.15, b i = 0.85 for the hippocampal surfaces. Correspondences between the shapes were established with the direct mapping techniques of [START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF][START_REF] Shi | Direct mapping of hippocampal surfaces with intrinsic shape context[END_REF] and the shapes were parameterized with the methods of [START_REF] Praun | Spherical parametrization and remeshing[END_REF]. Figure 2 further illustrates the method with the calculation of the sample Fréchet mean shape of 8 hippocampal surfaces.

Energy and Localization

A pre-shape geodesic p(t) = (r t , v t ), 0 t 1, has parallel velocity field. In particular, it is traversed with constant speed ω, where ω is the length of p. Thus, the energy of p is E(p) = 1 0 ∂ t p(t) 2 p(t) dt = ω 2 . On the other hand, we can write the energy as

E(p) = m i=1 1 0 a i |∂ t r t i | 2 e ri + b i ∂ t v t i 2 e ri dt , (12)
which expresses E(p) as a sum of the local contributions of the edges. Thus, the energy density function

ρ(e i ) = 1 ω 2 1 0 a i |∂ t r t i | 2 e ri + b i ∂ t v t i 2 e ri dt (13)
quantifies the fraction of the total deformation energy associated with the ith edge. Although geodesic distance is a global quantifier of shape difference, the energy density provides a means to measure local shape differences and identify the regions where shape similarity and divergence are most pronounced. One can further decompose the local energy into its tension and rigidity components to separately quantify the local shape differences due to stretching and bending.

One may also modify ρ to a function defined on the vertex set by letting the value on a vertex be the average value of ρ on the edges incident with that vertex. Figures 3 displays the initial shapes of the geodesics in Figure 1 overlaid with color maps of the respective energy density functions.

Summary and Discussion

We constructed a shape space of elastic spherical surfaces immersed in Euclidean space equipped with a family of geodesic metrics that depend on the elasticity of the surfaces. The metrics reflect the resistance offered by the shapes to deformations by stretching and bending. Although the general model is anisotropic, a special homogeneous sub-collection associated with constant tension and rigidity coefficients should be of particular interest in applications. The selection of elasticity parameters for a particular problem based on shape discrimination or other criteria requires further study. We developed an algorithm to calculate geodesics between spherical shapes and introduced energy density functions that provide a localization tool that allows us to identify the regions where shape dissimilarity is most pronounced. Energy density functions also let us characterize the nature of the local elastic deformations as due to stretching or bending. We illustrated the applicability of the methods with the calculation of geodesic deformations between pairs of cortical and hippocampal surfaces segmented from MR images. Geodesic distances, geodesic interpolations, and models of shape variation developed in shape space enable us to quantify and visualize anatomical resemblance and divergence between individuals and across populations. We applied the methods developed in the paper to the construction of an atlas of the hippocampus as the sample Fréchet mean shape of a group of 8 segmented hippocampi. Further statistical modeling of the shape of surfaces, extensions of the model to non-spherical shapes, and applications of the computational models to problems in neuroimaging are some of the topics to be investigated in future work.

A Closest-Shape Projection

Given (r, v) ∈ M near the pre-shape space P , we need an efficient estimation of the pre-shape closest to (r, v). Consider the kn + 1 residual functions ρ 1 (r, v) = 1 -G 1 (r, v) and ρ ℓ,j (r, v) = -G ℓ,j (r, v), whose vanishing is equivalent to (r, v) ∈ P . We estimate the closest pre-shape using Newton's method to find the zeros of H

(r, v) = 1 2 ρ 2 1 (r, v) + 1 2 ℓ,j ρ 2 ℓ,j (r, v).
If we relax the condition that each row of v must satisfy v i = 1 and just treat each v i as an arbitrary vector in R k , then the negative gradient of H is

-∇H(r, v) = ρ 1 (r, v)∇G 1 (r, v) + ℓ,j ρ ℓ,j (r, v)∇G ℓ,j (r, v) , (14)
where ∇G r 1 (r, v) = [1/a 1 . . . 1/a m ] T is the modular component of ∇G 1 (r, v) and the directional component is zero. The modular part of ∇G ℓ,j (r, v) and the jth column of its directional component, which is its only nonzero column, are 0

ε ℓ 1 a ℓ 1 v ℓ1,j 0 ε ℓ 2 a ℓ 2 v ℓ2,j 0 ε ℓ 3 a ℓ 3 v ℓ3,j 0 T and 0 ε ℓ 1 b ℓ 1 0 ε ℓ 2 b ℓ 2 0 ε ℓ 3 b ℓ 3 0 T . (15)
Here, the nonzero entries occur in rows ℓ 1 , ℓ 2 and ℓ 3 . Write the modular and directional parts of the gradient as -∇H(r, v) = (δ r , δv ). To make δv i tangential to S 2 at v i , replace it with

δ v i = δv i -( δv i • v i )v i . Letting ǫ(r, v) = H(r, v)/ (δ 1 , δ 2 ) 2
(r,v) , update (r, v) as follows:

r = r + ǫδ 1 ; v i = cos(ǫ δ v i )v i + sin(ǫ δ v i ) δ v i δ v i , if δ v i = 0 ; (16)
and v i stays unchanged, otherwise. Iterate until H(r, v) becomes small.

B Covariant Integration

For (z, y) = (z, y 1 , . . . , y k ) ∈ R × R k , we use the subscript 0 to identify the z-coordinate and the subscript j, 1 j k, for the coordinate y j . Given a i , b i > 0, define a Riemannian metric on R × R k whose metric tensor at (z, y) is g 00 (z, y) = a i e z , g jj (z, y) = b i e z , and 0, otherwise. Then, the Riemannian structure on E, defined in Section 3, is isometric to the Cartesian product of these (k + 1)-dimensional models over 1 i m. Thus, to derive the differential equation that governs covariant integration along a path in M , it suffices to derive the corresponding differential equation for covariant integration in R × S k-1 with respect to the induced metric. The Christoffel symbols of the Levi-Civita connection on R × R k are Γ 0 00 = Γ j 0j = Γ j j0 = 1/2, Γ 0 jj = -b i /2a i , 1 j k, and zero otherwise. Therefore (cf. [START_REF] Do Carmo | Riemannian Geometry[END_REF]), the covariant derivative of a vector field (F t , X t ) along a path (z t , y t ) in R × R k is given by

D t F t = ∂ t F t + 1 2 (∂ t z t )F t -1 2 bi ai [X t • ∂ t y t ] D t X t = ∂ t X t + 1 2 (X t ∂ t z t + F t ∂ t y t ) . (17)
If (z t , y t ) is a path in R × S k-1 , then fields that are tangential to R × S k-1 are those that satisfy the additional orthogonality condition F t • y t = 0. Thus, if (f t , x t ) and (F t , X t ) are both tangential to R × S k-1 , the fact that the covariant derivative of (F t , X t ) in the submanifold R × S k-1 is (f t , x t ) may be rephrased as D t F t = f t and D t X t = x t + τ t y t , where τ t is a scalar field to be determined. Substituting in (17), we obtain

∂ t F t = f t -1 2 ∂ t z t F t + 1 2 bi ai [X t • ∂ t y t ] ∂ t X t = x t -1
2 (X t ∂ t z t + F t ∂ t y t ) + τ t y t .

(18) Differentiating X t • y t = 0, we get ∂ t X t • y t = -X t • ∂ t y t . From (18), it follows that τ t = -X t • ∂ t y t , where we used the facts that X t • y t = 0 and ∂ t y t • y t = 0. Substituting this value of τ t in (18), we obtain [START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF].

Fig. 1 .

 1 Fig. 1. (a) Two cortical surfaces extracted from MR images and three intermediate stages of a geodesic interpolation with homogeneous coefficients ai = 0.3 and bi = 0.7; (b) a similar calculation with hippocampal surfaces with ai = 0.15 and bi = 0.85.

Fig. 2 .

 2 Fig. 2. The mean of 8 right hippocampal surfaces calculated with homogeneous coefficients ai = 0.15 and bi = 0.85. The mean shape is shown on the last panel.

Fig. 3 .

 3 Fig. 3. Plots of the energy density functions for the geodesics shown in Figure 1.

⋆ This work was supported in part by NSF grants DMS-0713012 and CCF-0514743, and NIH Roadmap for Medical Research grant U54 RR021813.


Grace Vesom

Nathan D Cahill

Lena Gorelick

J Alison Noble

Characterization of Anatomical Shape Based on Random Walk Hitting Times

This paper presents an implicit shape representation for describing anatomical shapes with high inter-patient variability based on the expected boundary hitting time of a random walk, which happens to be the solution to the Poisson equation. The main contribution of this paper is to test the validity of the Poisson-based mapping for learning anatomical shape variability, comparing its compactness and completeness with the commonly used Signed Distance Transform and using the liver and the caudate nucleus as examples. Based on these findings, we discuss its use as a shape prior for image segmentation.

Introduction

An integral part of modern disease management is treatment planning, which involves several aspects of medical image analysis, from segmenting objects of interest for volume measurement to intra-patient registration for monitoring size change or morphology. Because of the difficulty of these tasks, many segmentation or registration frameworks require a shape representation for the object in question. Shape representations can be abstractly divided into two categories: explicit and implicit. Explicit shape representations are parameterized to create ordered collections of components. Examples include Active Shape Models [START_REF] Cootes | Active shape models-their training and application[END_REF] and spherical harmonics [START_REF] Brechbuhler | Parametrization of closed surfaces for 3-D shape description[END_REF]. Mostly, they define a (N -1)-dimensional boundary in a N -dimensional image space, but can include additional information such as a medial representation [START_REF] Styner | Medial models incorporating object variability for 3D shape analysis[END_REF] or, like Active Appearance Models [START_REF] Cootes | Active appearance models[END_REF], intensity and texture.

An implicit representation of an N -dimensional shape is a function f (x) in an N -dimensional space that takes on the value zero if and only if x is a position on the (N -1)-dimensional shape boundary. The function is defined over the inside domain of the closed object, and sometimes, depending on the application, over the entire image domain.

Implicit shape descriptions are free of topological constraints and have local support. In recent literature, the Signed Distance Transform (SDT) has proved to be a popular implicit shape representation -attractive for ease of computation, definition over the entire image space, and ability to fit into a level-set segmentation framework. In addition to the SDT, Hong et al. [START_REF] Hong | Shape Representation based on Integral Kernels: Application to Image Matching and Segmentation[END_REF] developed an implicit shape representation based on an approximation to the solution of the heat equation over the shape. Their shape representation is attractive because it can be quickly computed by convolution with a Gaussian kernel. However, it does not adjust the Gaussian kernel to account for arbitrary shape boundaries.

To be able to summarize anatomical object variability, we seek a shape descriptor which can extract characteristic global object structure while preserving individual variations. In our work, we consider another implicit shape representation called the Poisson Transform (PT) and look at its application in learning and modeling anatomical variation. The PT was first used by Gorelick et al. [START_REF] Gorelick | Shape Representation and Classification Using the Poisson Equation[END_REF] to analyze 2D shape properties on silhouettes in order to perform shape classification and retrieval. We extended that idea in [START_REF] Cahill | Investigating Implicit Shape Representations for Alignment of Livers from Serial CT Examinations[END_REF] to include the shape exterior and used this in the context of deformable image registration. In the next section, we revisit the mathematical representation of the SDT, then provide a summary of and computational method for the PT. We use this new mapping to produce an implicit shape representation, and then analyze it for completeness and compactness in comparison to the SDT. We use Principal Component Analysis on the SDT for our training data to summarize shape variation in the same manner as [START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF][START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF][START_REF] Rousson | Implicit Active Shape Models for 3D Segmentation in MR Imaging[END_REF]. We investigate the validity of the SDT representation for anatomical objects with highly variable shape against the PT representation. Besides the technical report by Lamecker et al. [START_REF] Lamecker | Segmentation of the Liver using a 3D Statistical Shape Model[END_REF], no studies have assessed the compactness or completeness of liver shape description. We conclude by discussing its application to image segmentation.

Implicit Shape Representations

Let us define S to be the interior of the shape and assume that S ⊂ Ω, where Ω is the entire image space. Then ∂S = S -S is the boundary of S, where S represents the closure of the open set S.

Signed Distance Transform

The SDT, also known as the signed Euclidean distance transform, is represented here by W SDT . It yields two pieces of information; the magnitude provides the Euclidean distance between a point and the closest point on ∂S, and the sign indicates whether the current point is inside S (negative) or outside S (positive). Mathematically speaking, W SDT is a solution to the Eikonal equation:

|∇W SDT (x)| = 1 ∀x ∈ Ω, W SDT (x) = 0 ∀x ∈ ∂S. (1)
Examples of the SDT for the liver and caudate nucleus can be seen in Figs. 1(a) and 1(c).

Random Walk Hitting Time

Gorelick et al. [START_REF] Gorelick | Shape Representation and Classification Using the Poisson Equation[END_REF] presented an implicit shape representation based on the expected time for a symmetric random walk to reach the shape boundary. Based on the solution to Poisson's equation, this representation has several advantages as it is differentiable everywhere and has the interior regularity property.

In order to derive the Poisson Transform, consider a shape S defined on a 3-dimensional discrete lattice with isotropic spacing h. Now consider a particle initially at position (x, y, z) ∈ S. If the particle undergoes a symmetric random walk, then we define U (x, y, z) to be the expected time (in number of steps) for the particle to reach any point on the boundary of S, given that the particle started at (x, y, z). This expected time is also known as the "hitting time" because it represents the average time it will take the particle to "hit" the boundary. Therefore, if (x, y, z) ∈ ∂S, then U (x, y, z) = 0. If (x, y, z) ∈ S, then the hitting time at (x, y, z) can be related to the hitting time at each of its six neighbors by a conditioning argument:

U (x, y, z) = 1 + 1 6 U (x + h, y, z) + U (x -h, y, z) U (x, y + h, z) + U (x, y -h, z) (2)
U (x, y, z + h) + U (x, y, z -h) .
Note that ( 2) is a discretized version of the Poisson equation

∆U (x, y, z) = - 6 h 2 , (3)
with ∆U = U xx + U yy + U zz denoting the Laplacian of U . For simplicity, we chose h = 1.

To use comparable notation with the SDT definition, let W PT replace U to represent the symmetric random walk hitting time to the shape boundary in an open domain, and let x = (x, y, z) for the remainder of this paper. Like the SDT, we seek to use the sign of the metric values to define points inside and outside the boundary shape, which results in changing the sign of (3), while the magnitude of those values still provides the symmetric random walk time to a point on the boundary. Thus, for the shape interior, the Poisson Transform W PT satisfies

∆W PT (x) = 6 ∀x ∈ S, W PT (x) = 0 ∀x ∈ ∂S. (4)
On the exterior of the shape, if no further boundary conditions are prescribed, the random walk analogy fails and the Poisson equation has infinitely many solutions. Therefore, in order to provide an extension of this shape representation to the entirety of Ω, we need to define some sort of external boundary condition. One option is to enforce Neumann boundary conditions on ∂Ω; however, this choice gives undesirable behavior of the gradient of the Poisson Transform at the image boundaries. Instead, we choose to define an open sphere T ⊃ Ω centered on the centroid of Ω, and solve the following boundary value problem:

∆W PT (x) = -6 ∀x ∈ T -S, ∇W PT (x) , n(x) = 0 ∀x ∈ ∂T, (5)
W PT (x) = 0 ∀x ∈ ∂S,
where n(x) is the outward pointing normal vector to the surface ∂T at x. The Neumann boundary condition on ∂T effectively reflects the symmetric random walk in the direction normal to the sphere. Once W PT has been found on T , it can be cropped to Ω.

Computing a discrete approximation of W PT inside S can be done in linear time with a multigrid algorithm [START_REF] Gorelick | Shape Representation and Classification Using the Poisson Equation[END_REF]. This algorithm can be extended to approximate W PT inside T -S by appropriately discretizing the Neumann boundary conditions, as established by Cahill et al. [START_REF] Cahill | Investigating Implicit Shape Representations for Alignment of Livers from Serial CT Examinations[END_REF]. We have found that applying two w-cycles that employ naive boundary conditions at the coarse levels (placing the boundary at the nearest coarse grid points, instead of modifying the nearby coarse equations to account for the fine, pixel-level location of the boundary) provides a good approximation. For further details on the multigrid method, see [START_REF] Trottenberg | Multigrid[END_REF].

Examples of the SDT for the liver and caudate nucleus can be seen in Figs.

Application to Segmentation

In medical image segmentation, image data information is typically balanced by integrating a priori knowledge into the algorithm through a shape prior, providing a global constraint in a segmentation framework.

Shape Representation

Using a shape prior for image segmentation requires an application-specific choice of shape representation. Selecting a shape representation or descriptor is a function of the information to be extracted from an image, the signal-tonoise of the image, the anatomical object in question, and the ultimate goal of the task. For the purpose of segmentation, the shape descriptor also subsequently affects the dissimilarity measure for comparing shapes and method of integration for a model into a segmentation framework. Golland et al. [START_REF] Golland | Detection and analysis of statistical differences in anatomical shape[END_REF] found that the choice of shape representation changed their classifier function results by about 15%, independent of shape analysis and alignment.

Shape Analysis Methods

With highly variable anatomical shapes, large numbers of samples, or training shapes, are required to correctly portray a shape population. Shape models are used to efficiently summarize the training shape space. Much research looks for the best shape analysis method to summarize the shape space, while some methods never question the validity of its use, such as with Principal Component Analysis (PCA). That method is utilized to reduce the dimensionality of the training shape space, while preserving the most variation. PCA is often challenged because of its sensitivity to low sample size, as well as the assumption of a Gaussian uni-modal distribution for the features measured in the resulting model. In the case of the SDT and PT, the shape representation space is not closed under linear operations, but PCA is still often employed. Shape analysis alternatives include independent component analysis [START_REF] Hyvarinen | Independent component analysis: algorithms and applications[END_REF], principal factor analysis [START_REF] Ballester | On the adequacy of principal factor analysis for the study of shape variability[END_REF], and support vector machines [START_REF] Golland | Detection and analysis of statistical differences in anatomical shape[END_REF]. Additionally, methods to work around obstacles in PCA have shown some success. They transform training data into a feature space via Mercer kernels [START_REF] Cremers | Shape statistics in kernel space for variational image segmentation[END_REF], or into a vector space using Fourier approximation [START_REF] Srivastava | Statistical shape analysis: clustering, learning, and testing. Pattern Analysis and Machine Intelligence[END_REF] or the logarithm of odds ratio [START_REF] Pohl | Using the logarithm of odds to define a vector space on probabilistic atlases[END_REF], where linear operations of shape representation transformations are closed.

With the SDT, PCA has been consistently employed for shape analysis for the purpose of learning shape variability in creating a statistical shape model. Notwithstanding, Hong et al. [START_REF] Hong | Shape Representation based on Integral Kernels: Application to Image Matching and Segmentation[END_REF] used their integral kernel shape representation as a template with a deformation field for segmentation, removing the need for any shape analysis.

Implicit Anatomical Shape Priors

In image segmentation of highly variable shapes, it is essential to integrate a priori information as a global constraint with data-driven energy functions. We concentrate here on statistical shape models (SSMs), as opposed to deformable shape models or classifiers.

Leventon et al. [START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF] employed SSMs from implicit shape representations, using the SDT to represent shapes and PCA to summarize the shape space to create a shape prior in image segmentation. Thereafter, literature involving implicit shape models concentrates on shape alignment prior to PCA and shape prior integration into a region-based level-set segmentation framework. Rousson et al. [START_REF] Rousson | Implicit Active Shape Models for 3D Segmentation in MR Imaging[END_REF] used the sum of squared distances as a dissimilarity measure between the prior and target, while minimizing the rotation and translation pose parameters during the level set evolution. Tsai et al. [START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF] similarly proposed an alignment method, but also included scale as a pose parameter and explored three energy functionals, optimizing the functionals via gradient descent. Yang et al. [START_REF] Yang | Neighbor-constrained segmentation with level set based 3-D deformable models[END_REF] took another step by introducing neighbor priors, developing a hierarchical multiobject segmentation. Accounting for neighboring structures in the shape prior, however, gave rise to concerns about closedness of linear combination of signed distance functions. Neighboring structures were modeled by the difference of their representative SDTs, however any linear combination of these differences may not represent neighboring structures as mutually exclusive.

Analysis of Implicit Shape Representations SDT and PT

While any dimensionality reduction method can arguably be more suitable for shape analysis, we wish to show that a shape descriptor that suppresses redundancy in variability, yet elicits distinguishing global features, will result in more meaningful factors. PCA is the most widely used form of shape analysis in medical imaging because of its ease of computation and interpretation. For these reasons, we have thus chosen to use PCA for a comparative shape analysis study on the SDT and PT.

The nature of the PT allows us to combine global and local shape extraction -we can see a common shape for each structure, while local variations are preserved closer to the boundary (Figs. 3(b) and 2(b)). Consequently, we believe PT is best suited to represent objects with high shape variability, for which the liver and caudate nucleus are well known and well studied. Consider the case of representing the shape of the liver. Although liver shape has consistent convexities, concavities, and ridges, it is still very highly variable across the population. It is also adjacent to structures of similar density, making it extremely difficult to segment in Computed Tomography (CT) images. Liver image segmentation has been a well-established problem in medical imaging analysis. A variety of published methods range from deformable models using the evolutionary algorithm [START_REF] Heimann | A Shape-Guided Deformable Model with Evolutionary Algorithm Initialization for 3D Soft Tissue Segmentation[END_REF] and data-driven methods [START_REF] Soler | Fully automatic anatomical, pathological, and functional segmentation from CT scans for hepatic surgery[END_REF] to graph-cuts [START_REF] Boykov | Graph Cuts and Efficient N-D Image Segmentation[END_REF], demonstrating the difficulty of the task.

The shape of the caudate nucleus is often investigated in order to identify brain abnormalities. An irregular shape in longitudinal studies and intra-patient asymmetric shape can reveal developmental disorders, such as schizophrenia [START_REF] Levitt | Shape of caudate nucleus and its cognitive correlates in neuroleptic-naive schizotypal personality disorder[END_REF].

In 3D, we computed the SDT and PT for 33 liver images and 20 caudate nucleus images. Slices through the mapping for each transform and class are shown in Fig. 1, and 3D renderings of the inner level sets are shown in Figs. 3 and2.

The Data

Our data has been obtained from the MICCAI 3D Segmentation Challenge [START_REF] Van Ginneken | 3D Segmentation in the Clinic: A Grand Challenge[END_REF] and the Churchill Hospital, part of the Oxford Radcliffe Hospitals NHS Trust in Oxford, UK. We used the segmented training data sets from MICCAI, and we manually segmented the additional liver data sets subject to the approval of a clinician. All binary images were rescaled to physical Cartesian coordinates and smoothed to remove aliasing effects. The liver surfaces were aligned as binary images using the Principal Axis Transformation [START_REF] Alpert | The Principal Axes Transformation-A Method for Image Registration[END_REF], while the caudate nucleus shapes were aligned using first-order moments, with a few slight manual rotations.

Analysis

For each case, we ran PCA over the entire data set, Figs. 4(a We defined compactness of a shape description as requiring as few parameters as possible and show relative compactness here in comparison to the SDT. From Figs. 4(a) and 5(a), we see that to summarize 95% of shape variance for the liver and the caudate nucleus, the SDT requires 15 and 8 modes, while the PT requires only 5 and 3, respectively -approximately a 3-fold improvement in both cases.

Discussion

With a large number of training data, we would expect that the number of modes required should approach a constant value, as we would expect in a complete shape model. If the number of modes required grows with the number of training data, then it can be inferred that every randomly chosen set of 5 training data added were somehow orthogonal to the existing set, which is improbable. Then PCA might be incapable of accurately reducing the dimensionality of a shape space spanned by the SDT on liver and caudate nucleus shapes. As another possibility, the SDT might be an inadequate representation of highly variable For a SSM based on the liver for image segmentation, Lamecker et al. [START_REF] Lamecker | Segmentation of the Liver using a 3D Statistical Shape Model[END_REF] used a subsampled geometric mesh representation. Mesh point values, like all values in a SDT-map, are directly affected by any noise on the shape surface. Lamecker report compactness and completeness curves consistent with that of the signed distance function in their liver SSM analysis, and required 18 modes to summarize 95% variance for 33 training shapes and 21 modes for 42 training shapes.

In contrast, the PT completeness curve flattens quickly as the number of training shapes grows, converging to a fixed number of modes, as we should expect in a complete model (Fig. 4(b)). Through our comparative PCA study we can see the PT is a relatively more compact and complete shape descriptor than SDT.

For high dimensional data and low sample size data sets, PCA is argued to be incomplete. That is to say some actual shapes may be unrealizable from the basis, which is chosen for low-approximation while summarizing the most variance. Yet one can argue that the shape analysis results are dependent on what has been extracted from the original data into the shape representation. Boundary representations give information about the object delineation, and SDT-maps propagate boundary variations beyond the neighborhood of occurrence, creating superfluous variations throughout the image domain. However, the PT-maps reveal global shape information from a data set yet adapt to individual object boundaries (see Figs. 12).

Conclusion and Future Work

In this paper, we introduced using a new implicit shape representation called the Poisson Transform, an extrapolation of Gorelick's innovative approach to shape characterization. In our comparative dimensionality reduction study, we showed that the PT gives approximately a 3-fold improvement over the SDT in the number of principal components required to represent 95% of the liver and caudate nucleus data set variance. The PT is capable of local support on the boundary for individual instances but substantially stable against those variations to give rise to a characteristic object shape across a data set. This balance in description enables optimized learning for highly variable anatomical objects. We are currently working on methods of integrating a shape model based on the PT into an image segmentation framework and overcoming obstacles mentioned in [START_REF] Golland | Detection and analysis of statistical differences in anatomical shape[END_REF][START_REF] Pohl | Using the logarithm of odds to define a vector space on probabilistic atlases[END_REF].

  1(b) and 1(d).

Fig. 1 .

 1 Fig. 1. Slices from example 3D liver and caudate nucleus images of both implicit functions, inside and outside the shape

Fig. 2 .Fig. 3 .

 23 Fig. 2. Inside level sets of the liver from SDT and PT

  ) and 5(a), and then for randomly chosen sets in multiples of 5, to show the number of modes needed to summarize a cumulative variance threshold of 95% given the changing number of training shapes, Figs. 4(b) and 5(b).

Fig. 4 .

 4 Fig. 4. Compactness and completeness of the SDT and PT for the shape of the liver

Fig. 5 .

 5 Fig. 5. Compactness and completeness of the SDT and PT for the shape of the caudate nucleus

Acknowledgments

The authors would like to thank Professor J. M. Brady for his insightful discussions; Tobias Heimann, Martin Styner, and Bram van Ginneken for allowing us to use their data; and Dr. Joanne Brady from the Churchill Hospital in Oxford for selecting patients and overseeing manual segmentations.


Moo K Chung
email: mkchung@wisc.edu

Anqi Qiu

M Brendon

Nacewicz

Seth Pollak

Richard J Davidson

Tiling Manifolds with Orthonormal Basis

One main obstacle in building a sophisticated parametric model along an arbitary anatomical manifold is the lack of an easily available orthonormal basis. Although there are at least two numerical techniques available for constructing an orhonormal basis such as the Laplacian eigenfunction approach and the Gram-Smidth orthogonalization, they are computationally not so trivial and costly. We present a relatively simpler method for constructing an orthonormal basis for an arbitrary anatomical manifold. On a unit sphere, a natural orthonormal basis is the spherical harmonics which can be easily computed. Assuming the manifold is topologically equivalent to the sphere, we can establish a smooth mapping ζ from the manifold to the sphere. Such mapping can be obtained from various surface flattening techniques. If we project the spherical harmonics to the manifold, they are no longer orthonormal. However, we claim that there exists an orthonormal basis that is the function of spherical harmonics and the spherical mapping ζ. The detailed step by step procedures for the construction is given along with the numerical validation using amygdala surfaces as an illustration. As an application, we propose the pullback representation that reconstructs surfaces using the orthonormal basis obtained from an average template. The pullback representation introduces less inter-subject variability and thus requires far less number of coefficients than the traditional spherical harmonic representation.

Introduction

We present a novel orthonormal basis construction method for an arbitrary anatomical surface that is topologically equivalent to a sphere. The method avoids the well known Gram-Smidth orthogonalization procedure [START_REF] Gorski | On determining the spectrum of primordial inhomogeneity from the cobe dmr sky maps: I. method[END_REF], which is inefficient for high resolution polygonal meshes. In order to perform the Gram-Smidth orthogonalization as described in [START_REF] Gorski | On determining the spectrum of primordial inhomogeneity from the cobe dmr sky maps: I. method[END_REF], for a surface mesh with n vertices, we need to perform the Choleski decomposition as well as the inversion of matrix of size n×n. For a cortical mesh generated with FreeSurfer [START_REF] Fischl | Measuring the thickness of the human cerebral cortex from magnetic resonance images[END_REF], n can easily reach up to 200000.

On the other hand, Qiu et al. [START_REF] Qiu | Smooth functional and structural maps on the neocortex via orthonormal bases of the laplace-beltrami operator[END_REF] constructed an orthonormal basis as the eigenfunctions of the Laplace-Beltrami operator in a bounded regions of interest (ROI) on a cortical surface (Figure 3). The finite element method (FEM) is used to numerically construct the orthonormal basis by solving a system of large linear equations. The weakness of the FEM approach is the computational burden of inverting a matrix of size n × n. One advantage of the eigenfunction approach is that since the eigenfunctions and eigenvalues are directly related to the Laplace-Beltrami operator, it is trivial and geometrically intuitive to construct the heat kernel analytically and perform a various heat kernel smoothing based modeling [START_REF] Chung | Cortical thickness analysis in autism with heat kernel smoothing[END_REF].

We propose a completely different method that avoids the computational bottleneck by using a conceptually different machinery. We assume an arbitrary anatomical surface to be topologically equivalent to a sphere. Then using a smooth mapping ζ obtained from a surface flattening technique, we project the spherical harmonics to the anatomical surface. Obviously the projected spherical harmonics will no longer be orthonormal. However, if we correct the metric distortion introduced from the surface flattening, we may able to make the projected spherical harmonics orthonormal somehow. This is the basic idea behind our new proposed method. For the surface flattening, we present a new method that treats the mapping ζ as the geodesic path of the heat equilibrium state.

As an application of the proposed technique, we present the novel pullback representation for parameterizing anatomical boundaries that outperforms the traditional spherical harmonic (SPHARM) representation [START_REF] Brechbuhler | Parametrization of closed surfaces for 3d shape description[END_REF] [3] [8] [START_REF] Shen | surface-based approach for classification of 3d neuroanatomical structures[END_REF] [START_REF] Styner | Framework for the statistical shape analysis of brain structures using spharm-pdm[END_REF]. We claim that our proposed representation has far less intersubject variability in the estimated parameters than SPHARM and converges faster to the true boundary with less number of basis.

Methods

It is assumed that the anatomical boundary M is a smooth 2-dimensional Riemannian manifold parameterized by two parameters. The one-to-one mapping ζ from point p = (p 1 , p 2 , p 3 ) ′ ∈ M to u = (u 1 , u 2 , u 3 ) ′ ∈ S 2 , a unit sphere, can be obtained from various surface flattening techniques such as conformal mapping [START_REF] Angenent | On the laplace-beltrami operator and brain surface flattening[END_REF] [8] [START_REF] Hurdal | Cortical cartography using the discrete conformal approach of circle packings[END_REF], quasi-isometric mapping [START_REF] Timsari | An optimization method for creating semi-isometric flat maps of the cerebral cortex[END_REF], area preserving mapping [START_REF] Brechbuhler | Parametrization of closed surfaces for 3d shape description[END_REF] [16] [START_REF] Styner | Framework for the statistical shape analysis of brain structures using spharm-pdm[END_REF] and the deformable surface algorithm [START_REF] Macdonald | Automated 3-D extraction of inner and outer surfaces of cerebral cortex from mri[END_REF]. Since the conformal mapping tend to introduce huge area distortion, most spherical harmonic literature tend to use area preserving mapping [START_REF] Brechbuhler | Parametrization of closed surfaces for 3d shape description[END_REF] [16] [START_REF] Styner | Framework for the statistical shape analysis of brain structures using spharm-pdm[END_REF].

In this paper, we present a new flattening technique via the geodesic trajectory of the equilibrium sate of heat diffusion. The proposed flattening technique is numerically simpler than any other available methods and does not require optimizing a cost function. The methodology is illustrated using the 47 amyg-Fig. 1. The diffusion equation with a heat source (amygdala) and a heat sink (enclosing sphere) corresponds. After sufficient amount of diffusion, the heat equilibrium state is reached. By tracing the geodesic path from the heat source to the heat sink using the geodesic contours, we obtain a smooth mapping ζ. dala binary segmentation obtained from the 3-Tesla magnetic resonance images (MRI).

High resolution anatomical MRI were obtained using a 3-Tesla GE SIGNA scanner with a quadrature head coil. Details on image acquisition parameters are given in [START_REF] Nacewicz | Amygdala volume and nonverbal social impairment in adolescent and adult males with autism[END_REF]. MRIs are reoriented to the pathological plane for optimal segmentation and comparison with an atlas. This global alignment guarantee that amygdala are approximately aligned in the same orientation.

Manual amygdala segmentation was done by a trained expert and the reliability of the manual segmentation was validated by two raters on 10 amygdale resulting in interclass correlation of 0.95 and the intersection over the union of 0.84 [START_REF] Nacewicz | Amygdala volume and nonverbal social impairment in adolescent and adult males with autism[END_REF]. Afterwards a marching cubes algorithm was used to extract the boundary of the binary segmentation as a triangle mesh with approximately 2000-3000 vertices. The amygdala surface is then mapped onto a sphere using the new flattening algorithm.

Diffusion-Based Surface Flattening

Given an amygdala binary segmentation M a , we put a larger sphere M s that encloses the amygala (Figure 1 left). The amygdala and sphere serve as Dirichlet boundary conditions for solving the Laplace equation. The amygdala is assigned the value 1 while the enclosing sphere is assigned the value -1, i.e.

f (M a , σ) = 1, f (M s , σ) = -1 (1)
for all σ. The amygdala and the sphere serve as a heat source and a heat sink respectively. Then we solve an isotropic diffusion

∂f ∂σ = ∆f (2)
within the empty space bounded by the amygdala and the sphere. ∆ is the 3D Laplacian. After enough diffusion, the system reaches the heat equilibrium state where the additional diffusion does not make any difference in the heat distribution (Figure 1 middle). Once we obtained the equilibrium state, we trace the geodesic path from the heat source to the heat sink for every mesh vertices. The trajectory of the geodesic path provides a smooth mapping from the amygdala surface to the sphere. The geodesic path can be easily traced by constructing geodesic contours that correspond to the level set of the equilibrium state (Figure 1 right). Then the geodesic path is constructed by finding the shortest distance from one contour to the next and iteratively connecting the path together. Figure 2 shows the process of flattening using five contours corresponding to the temperature 0.6, 0.2, -0.2, -0.6, -1.0. Although we did not apply our flattening technique to other anatomical objects, the proposed method can be applied to more complex object than the amygdala. At the equilibrium state, we no longer has change in heat change over time, i.e. ∂f ∂σ = 0, so we have the Laplace equation ∆f = 0 with the same boundary condition. The Laplace equation has been previously used to trace the distance between outer and inner cortical surfaces and to compute cortical thickness [START_REF] Jones | Three-dimensional mapping of cortical thickness using laplace's equation[END_REF] [12] [START_REF] Yezzi | An eulerian pde approach for computing tissue thickness[END_REF]. Since the solution to the Laplace equation with the boundary condition ( 1) is unique even for highly convoluted and folded structures, the geodesic path will be uniquely defined.

Orhonormal basis in two sphere S 2

Suppose a unit sphere S 2 is represented as a high resolution triangle mesh consisting of the vertex set V(S 2 ). We have used an almost uniformly sampled mesh with 2562 vertices and 5120 faces. Let us parameterize coordinates u ∈ S 2 with parameters θ, ϕ:

(u 1 , u 2 , u 3 ) = (sin θ cos ϕ, sin θ sin ϕ, cos θ),
where (θ, ϕ) ∈ N = [0, π] ⊗ [0, 2π). The polar angle θ is the angle from the north pole and ϕ is the azimuthal angle. The orthonormal basis on the unit sphere is given by the eigenfunctions of

∆f + λf = 0,
where ∆ is the spherical Laplacian. The eigenfunction Y lm corresponding to the eigenvalue l(l + 1) is called the spherical harmonic of degree l and order m [START_REF] Courant | Methods of Mathematical Physics[END_REF].

With respect to the inner product

f, g S 2 = S 2 f (u)g(u) dµ(u), (3)
with measure dµ(u) = sin θdθdϕ, Y lm form the orthonormal basis in L 2 (S 2 ), the space of square integrable functions on S 2 , i.e.

Y lm , Y l ′ m ′ S 2 = δ ll ′ δ mm ′ . (4)
The inner product can be numerically computed as the Riemann sum over mesh vertices as

Y lm , Y l ′ m ′ S 2 ≈ uj ∈V(S 2 ) Y lm (u j )Y l ′ m ′ (u j )D S 2 (u j ), (5)
where D S 2 (u j ) is the discrete approximation of dµ(u). Let T 1 uj , T 2 uj , • • • , T jm uj be the area of triangles containing the vertex u j . Then we estimate D S 2 (u j ) as

D S 2 (u j ) = 1 3 jm k=1 T k uj . (6)
The discrete approximation (6) defines the area of triangles at a mesh vertex. The factor 1/3 is chosen in such a way that uj ∈V(S 2 ) D S 2 (u j ) = 12.5514 = 4 • 3.1378, analogous to the relationship

S 2 dµ(p) = 4π.
The discrepancy between the integral and its discrete counter part is due to the mesh resolution and it should become smaller as the mesh resolution increases.

Based on the proposed discretization scheme, we have computed the inner product (5) for all degrees 0 ≤ l, l ′ ≤ 20. Figure 3 (left) shows the inner products for every possible pairs. Since for up to the k-th degree, there are total (k + 1) 2 basis functions, we have total 441 2 possible inner product pairs, which is displayed as a matrix. For the diagonal terms, we obtained 0.9988 ± 0.0017 while for the off-diagonal terms, we have obtained 0.0000 ± 0.0005 indicating our basis and the discretization scheme is orthonormal with two decimal accuracy.

Orthonormal basis on manifold M

For f, g ∈ L 2 (M), the orthonormality is defined with respect to the inner product

f, g M = M f (p)g(p) dµ(p).
Using the spherical harmonics in S 2 , it is possible to construct an orthonormal basis in M numerically without the computational burden of solving the large matrix inversion associated with the eigenfunction method or the Gram-Smidth orthogonalization. Since the spherical harmonics are orthonormal in S 2 and, the manifolds S 2 and M can be deformed to each other by the mapping ζ, one would guess that the orthonormal basis in M can be obtained somehow using the spherical harmonics. Surprisingly this guess is not wrong as we will show in this section.

For f ∈ L 2 (S 2 ), let us define the pullback operation * as

ζ * f = f • ζ.
While f is defined on S 2 , the pullbacked function ζ * f is defined on M. The schematic of the pull back operation is given in Figure 6 (a). Then even though we do not have orthonormality on the pullbacked spherical harmonics, i.e.,

ζ * Y lm , ζ * Y l ′ m ′ M = δ ll ′ δ mm ′ ,
we can make them orthonormal by using the Jacobian determinant of the mapping ζ somehow. Consider the Jacobian J ζ of the mapping ζ : p ∈ M → u ∈ S 2 defined as

J ζ = ∂u(θ, ϕ) ∂p(θ, ϕ) .
For functions f, g ∈ L 2 (S 2 ), we have the following change of variable relationship:

f, g S 2 = M ζ * f (p)ζ * g(p)| det J ζ | dµ(p). (7)
Similarly we have the inverse relationship given as

ζ * f, ζ * g M = S 2 f (u)g(u)| det J ζ -1 | dµ(u). (8)
By letting f = Y lm and g = Y l ′ m ′ in (7), we obtain

δ ll ′ δ mm ′ = M ζ * Y lm ζ * Y l ′ m ′ | det J ζ | dµ(p) (9)
Equation ( 9) demonstrates that functions

Z lm = | det J ζ | 1/2 ζ * Y lm (10)
are orthonormal in M. We will refer l as degree and m as order of the basis function. Using the Riesz-Fischer theorem [START_REF] Kolmogorov | Introductory real analysis[END_REF], it is not hard to show that Z lm form a complete basis in L 2 (M).

Numerical Implementation

Although the expression (10) provides a nice analytical form for an orthonormal basis for an arbitrary manifold M, it is not practical. If one want to use the basis [START_REF] Jones | Three-dimensional mapping of cortical thickness using laplace's equation[END_REF], the Jacobian determinant needs to be numerically estimated somehow. We present a new discrete estimation technique for the surface Jacobian determinant that avoids estimating unstable spatial derivative estimation. The Jacobian determinant J ζ of the mapping ζ can be expressed in terms of the Riemannian metric tensors associated with the manifolds S 2 and M. Consider determinants det g S 2 and det g M of the Riemannian metric tensors associated with the parameterizations u(θ, ϕ) and p(θ, ϕ) respectively. Note that the integral of the area elements √ det g S 2 and √ det g M with respect to the parameter space N gives the total area of the manifolds, i.e.

N det g S 2 dµ(θ, ϕ) = 4π, N det g M dµ(θ, ϕ) = µ(M).
Then we have the relationship

| det J ζ -1 | = √ det g M √ det g S 2 , | det J ζ | = √ det g S 2 √ det g M .
Note that the Jacobian determinant det J ζ measures the amount of contraction or expansion in the mapping ζ from M to S 2 . So it is intuitive to have this quantity to be expressed as the ratio of the area elements. Consequently the discrete estimation of the Jacobian determinant at mesh vertex u j = ζ(p j ) is obtained as Then our orthonormal basis is given by

| det J ζ | ≈ D S 2 (u j ) D M (p j ) .
Z lm (p j ) = D S 2 (ζ(p j )) D M (p j ) ζ * Y lm (p j ). (11)
The numerical accuracy can be determined by computing the inner product

Z lm , Z l ′ m ′ M ≈ pj ∈V(M) Z lm (p j )Z lm (p j )D M (p j ). = pj ∈V(M) ζ * Y lm (p j )ζ * Y l ′ m ′ (p j )D S 2 (ζ(p j )) = uj ∈V(S 2 ) Y lm (u j )Y l ′ m ′ (u j )D S 2 (u j ) = Y lm , Y l ′ m ′ S 2
Since this is tautology, the order of the numerical accuracy in Z lm is identical to that of spherical harmonics given in the previous section. There is no need for additional validation other than given in the previous section. Hence we conclude that our basis is in fact orthonormal within two decimal accuracy. Figure 4 shows the result of our numerical procedure applied to the average amygdala surface template. The template surface is constructed by averaging the surface using the spherical harmonic correspondence given in [START_REF] Chung | Weighted fourier representation and its application to quantifying the amount of gray matter[END_REF].

We have also constructed an orthonormal basis on a cortical surface with more than 40000 mesh vertices (Figure 5). The diagonal elements in the inner product matrix are 0.9999 ± 0.0001 indicating that our basis is orthonormal within three decimal accuracy. As the mesh resolution increases, we expect to have increased accuracy. The proposed orthonormal basis construction methods avoids inverting matrix of size larger than 40000 × 40000 associated with the eigenfunction approach and the Gram-Smidth orthogonalization process.

Although the pattern of tiling in the eigenfunction approach (Figure 3) and the pullback based method (Figure 4) looks different, it can be shown that they are actually linearly dependent.

Application: Pullback Representation

As an application of the proposed orthonormal basis construction, we present a new variance reducing Fourier Series representation that outperforms the traditional spherical harmonic representation [START_REF] Brechbuhler | Parametrization of closed surfaces for 3d shape description[END_REF] [3] [8] [START_REF] Shen | surface-based approach for classification of 3d neuroanatomical structures[END_REF] [START_REF] Styner | Framework for the statistical shape analysis of brain structures using spharm-pdm[END_REF]. We will call this method as the pullback representation.

The spherical harmonic (SPHARM) representation models the surface coordinates with respect to a unit sphere as

p(θ, ϕ) = k l=0 l m=-l p 0 lm Y lm (θ, ϕ) (12)
where p 0 lm = p, Y lm S 2 are spherical harmonic coefficients, which can be viewed as random variables. The coefficients are estimated using the iterative residual fitting algorithm [START_REF] Chung | Weighted fourier representation and its application to quantifying the amount of gray matter[END_REF] that breaks a larger least squares problem into smaller ones in an iterative fashion. The MATLAB code for performing the iterative residual fitting algorithm for arbitrary surface mesh is given in http://www.stat.wisc.edu /∼mchung/softwares/weighted-SPHARM/weighted-SPHARM.hmtl. Note that all MRIs were reoriented to the pathological plane guaranteeing an approximate global alignment before the surface flattening to increase the robustness of the coefficient estimation. The shortcoming of the spherical harmonic representation is that the reconstruction is respect to a unit sphere that is not geometrically related to the original anatomical surface. On the other hand, the pullback representation will reconstruct the surface with respect to the average template surface reducing substantial amount of variability compared to SPHARM.

In the pullback representation, we represent the surface coordinates with respect to the template surface M as

p(θ, ϕ) = k l=0 l m=-l p 1 lm Z lm (θ, ϕ) (13)
with p 1 lm = p, Z lm M . Then we claim that the pullback representation has smaller variance in the estimated coefficients so that Var(p 1 lm ) ≤ Var(p 0 lm ).

The equality in ( 14) is obtained when the template M becomes the unit sphere, in which case the spherical mapping ζ collapses to the identity, and the inner products coincide. We have computed the sample standard deviation of Fourier coefficients for 47 subjects using the both representations. In average, the SPHARM contains 88% more intersubject variability compared to the pullback representation (Figure 6 right). This implies that SPHARM is an inefficient representation and requires more number of basis to represent surfaces compared to the pullback method.

Although the pullback method is more efficient, the both representations (12) and ( 13) converge to each other as k goes to infinity. We have computed the squared Euclidean distance between two representations numerically (Figure 7). In average, the difference is 0.0569 mm for 20 degree representation negligible for 1mm resolution MR. Figure 7 also visually demonstrate that the pullback representation converges to the true manifold faster than SPHARM again showing the inefficiency of the SPHARM representation.

Conclusion

We have introduced a computationally efficient way to construct an approximate orthonormal basis on an arbitrary manifold by pulling back the spherical harmonics to the manifold and accounting for the metric distortion using the Jacobian determinant. The proposed technique is very general so that it can be applicable to other types of anatomical manifolds. The constructed basis on an amygdala is used to show the new pullback representation that reconstruct the manifold as linear combination of the basis functions.

Fig. 2 .

 2 Fig. 2. Amygala surface flattening is done by tracing the geodesic path of the heat equilibrium state. The numbers corresponds to the different the geodesic contours. For simple shapes like amygde, 5 to 10 contours are sufficient for tracing the geodesic path.

Fig. 3 .

 3 Fig.3. Left: inner products of eigenfunctions of the Laplacian for every pairs[START_REF] Qiu | Smooth functional and structural maps on the neocortex via orthonormal bases of the laplace-beltrami operator[END_REF]. The pairs are rearranged from low to high degree. Right: representative eigenfunctions Ψj on the left amygdala template surface obtained by solving ∆Ψj + λjΨj = 0.

Fig. 4 .

 4 Fig.4. Left: inner products of spherical harmonics computed using formula (3) for every pairs. The pairs are rearranged from low to high degree and order. There are total (20 + 1) 2 = 441 possible pairs for up to degree 20. Right: representative orthonormal basis Z lm on the left amygdala template surface.

Fig. 5 .

 5 Fig. 5. Orthonormal basis Z lm on a cortical surface. The basis is projected on a sphere to show how the nonuniformity of the Jacobian determinant is effecting the spherical harmonics Y lm . The color scale is thresholded at ±0.003 for better visualization.

Fig. 6 .

 6 Fig. 6. Left: schematic showing how the pullback operation * is working. Point p ∈ M is mapped to u ∈ S 2 via our new flattening technique. As an illustration f = Y3,2 + 0.6Y2,1 is plotted on S 2 . The function f is pulled back onto M by ζ. Right: sample standard deviation of Fourier coefficients of for 47 subjects plotted over the index of basis. In average, the traditional SPHARM representation (black) has 88% more variability than the pull back method (red).

Fig. 7 .

 7 Fig. 7. Comparison of SPHARM and the pullback representations for degree 5 to 25. Red colored numbers are the average Euclidean distance between two representations in mm.