A Bayesian Generative Model for Surface Template Estimation

Abstract

3D surfaces are important geometric models for many objects of interest in image analysis and Computational Anatomy. In this paper, we describe a Bayesian inference scheme for estimating a template surface from a set of observed surface data. In order to achieve this, we use the geodesic shooting approach to construct a statistical model for the generation and the observations of random surfaces. We develop a mode approximation EM algorithm to infer the maximum a posteriori estimation of initial momentum 𝜇, which determines the template surface. Experimental results of caudate, thalamus, and hippocampus data are presented.

1. Introduction

3D surfaces are important geometric models for many objects of interest in image analysis and Computational Anatomy. For example, they are often used to represent the shape of 3D objects, the surface of human faces, and the boundaries of brain structures or of other human organs. Most data analysis methods in this domain are template-centered, and a proper estimation of a template plays an important role to obtain high quality results. This paper is devoted to the description of statistically supported template estimation method which is adapted to surface data sets.

Our approach will be to build a generative statistical shape model in which the template is a parameter. We will then estimate it using maximum likelihood. This model relies on the very natural setting for which an observed surface is a noisy version of a random deformation of the template. This is the most generic and most basic approach of the deformable template paradigm, even if we add a small refinement by including a prior distribution on the template, based on what we will call a hypertemplate. Even with this global scheme which is fairly simple, we will see that implementing it in the context of surfaces will constitute a significant theoretical and numerical challenge.

At the exception of the recent work of [1], this approach significantly differs from what has been mostly proposed in the literature, in which most of the methods compute templates as averages over specific common parametrizations of the surfaces (using, for example, the sphere as a parameter space [2]). Parametric representations, however, are limited by the fact that, because they are defined a priori and independently for each object, they cannot be assumed to suitably align important features in a given data set of surfaces (i.e., give similar coordinates to similar features in the surfaces). This usually results in oversmoothed template surfaces (which is the equivalent of getting blurry template images in the case of image averaging). In [1], a similar diffeomorphic transformation model is used, but, as we will see, our Bayesian construction will provide a well-specified template whereas [1] needs to rely to topologically unconstrained approximations to end up with a manageable template.

In addition to the references above, there have been several publications addressing the issue of shape averaging over a dataset, although most of them work with 3D volume data or landmark points set. In several cases, the average is based on metric properties of the space of shapes [3–7], and the template is computed as an intrinsic average, minimizing the sum of square distances to each element of the set (Fréchet mean). Such methods have been implemented in the context of diffeomorphic deformation models (which are also models of interest in the present paper) for landmark matching [8], for 3D average digital atlas construction [9], and to quantify variability in heart geometry [10]. Other definitions of the average, adapted to situations in which the data is corrupted by noise, have been proposed [11–13], based on variational approaches (but not relying on a generative statistical model).

Our approach to build a generative statistical shape model is reminiscent of the construction developed in [14] for linear models of deformations, and in [15] for large diffeomorphic deformations in 3D volume image averaging. Adapting these ideas to surfaces will however require new algorithms and numerical developments.

In order to present our model, we need to first provide some background materials and notation, describing in particular the geodesic shooting equations that we will use to generate deformed surfaces. We will then introduce a random statistical model describing the generation and the observations of random surfaces. We then develop a Mode Approximation EM algorithm for surface averaging, to estimate the template from observations. In the optimization part, we derive and implement a new variational scheme, which is also applicable to surface matching, providing an alternative approach to the one originally proposed in [16, 17]. Finally, we present and discuss experimental results on caudate, thalamus, and hippocampus data.

2. EPDiff for Surface Evolution

We will base our random shape model on the so-called EPDiff equation, which describes the evolution of deformable structures (like images, surfaces, or landmarks) under the action of groups of diffeomorphisms. It is a geodesic equation for a Riemannian metric on diffeomorphisms, and describes a momentum conservation law in the associated Hamiltonian system. The reader interested by the theory behind this equation can refer to [18–20], but most of this background will not be needed for the present paper, in which we will only use the specific form of the equations for surface evolution. The term EPDiff comes from its determination as an Euler-Poincaré equation in the group of diffeomorphisms, as introduced in [21]. One of its main interests here is that it provides a numerically stable, easily described, Hamiltonian evolution over diffeomorphisms, which will constitute an ideal support for our shape models.

The EPDiff equations describe the combined time evolution of a diffeomorphism, denoted 𝜙(𝑡,⋅) and of what can be interpreted as a momentum, denoted 𝑝(𝑡,⋅). The initial conditions are always 𝜙(0,𝑥)=𝑥 for 𝜙, and some initial value, 𝑝₀, for 𝑝. This initial momentum will be a key component of the statistical model that will be built later on.

Let us start with the simplest form of the equation, which assumes that 𝑝₀ is a vector-valued function over ℝ^𝑑, that is, 𝑝₀∶ℝ^𝑑→ℝ^𝑑. It involves a smoothing kernel, 𝐾, defined on ℝ³×ℝ³, a typical choice being𝐾(𝑥,𝑦)=cstexp(−‖𝑥−𝑦‖²2𝜏²). (1) Letting ∇₁𝐾 denote the gradient of 𝐾 with respect to its first variable, the corresponding EPDiff equation is𝑑𝜙𝑑𝑡(𝑡,𝑥)=_ℝ3𝐾(𝜙(𝑡,𝑥),𝜙(𝑡,𝑦))𝑝(𝑡,𝑦)𝑑𝑦,𝑑𝑝𝑑𝑡(𝑡,𝑥)=−_ℝ3∇₁𝐾(𝜙(𝑡,𝑥),𝜙(𝑡,𝑦))(𝑝(𝑡,𝑥)⋅𝑝(𝑡,𝑦))𝑑𝑦. (2) Here, the notation 𝑎⋅𝑏 refers to the usual dot product between vectors in ℝ³. If 𝐾 is smooth enough, this system has solutions for arbitrary large 𝑡, and the mapping 𝑥↦𝜙(𝑡,𝑥) is a diffeomorphism at all times.

The interesting fact about these equations is that they can have singular variants that are described in a similar way and have the same existence properties. The simplest way to relate the variants to the previous equation is to replace the Lebesgue's measure in the integrals by another, possibly singular, measure. For example, taking a surface 𝑆₀ in ℝ³, we can use the volume form on 𝑆₀ as a reference measure and obtain the equations (in which 𝑝₀ and 𝑝(𝑡,⋅) need only to be defined on 𝑆₀)𝑑𝜙𝑑𝑡(𝑡,𝑥)=_𝑆0𝐾(𝜙(𝑡,𝑥),𝜙(𝑡,𝑦))𝑝(𝑡,𝑦)𝑑𝜔_𝑆0(𝑦),𝑑𝑝𝑑𝑡(𝑡,𝑥)=−_𝑆0∇₁𝐾(𝜙(𝑡,𝑥),𝜙(𝑡,𝑦))×(𝑝(𝑡,𝑦)⋅𝑝(𝑡,𝑦))𝑑𝜔_𝑆0(𝑦), (3) where 𝑑𝜔_𝑆0 is the volume form on 𝑆₀. Note that the first equation is defined for 𝑥∈ℝ³, but it suffices to use it for 𝑥∈𝑆₀ to obtain an equation for the evolving surface 𝑆_𝑡=𝜙(𝑡,𝑆₀)={𝑦=𝜙(𝑡,𝑥),𝑥∈𝑆₀}. (4)

We can write a discrete form of the equations by replacing 𝑑𝑦 by a sum of Dirac measures (at points 𝑥₁,…,𝑥_𝐿 in ℝ³), which gives, letting 𝑎_𝑙(𝑡)=𝑝(𝑡,𝑥_𝑙),𝑑𝜙𝑑𝑡(𝑡,𝑥)=^𝐿_𝑙=1𝐾(𝜙(𝑡,𝑥),𝜙(𝑡,𝑥_𝑙))𝑎_𝑙,𝑑𝑎_𝑘𝑑𝑡=−^𝐿_𝑙=1(𝑎_𝑙⋅𝑎_𝑘)∇₁𝐾(𝜙(𝑡,𝑥_𝑘),𝜙(𝑡,𝑥_𝑙)). (5) Similarly to (3), the first equation is valid for all 𝑥∈ℝ³, but it suffices to solve it for 𝑥=𝑥_𝑙, 𝑙=1,…,𝐿 to obtain the evolution of the point set 𝑥_𝑙(𝑡)=𝜙(𝑡,𝑥_𝑙). (6) Also, (5) can be seen as a discretization of (3) in which 𝑥₁,…,𝑥_𝐿 are the vertices of a triangulation of 𝑆₀, and 𝑎_𝑙(𝑡)=𝑝(𝑡,𝑥_𝑙)𝛿𝜎_𝑆0(𝑥_𝑙), where 𝛿𝜎_𝑆0(𝑥_𝑙) is the area of a surface element around 𝑥_𝑙.

The evolution of point sets is the most important from a practical point of view, since it is an ODE that can be easily implemented. Assuming a radial kernel 𝐾(𝑥,𝑦)=𝛾(‖𝑥−𝑦‖²) like in (1), and denoting 𝛾_𝑘𝑙=𝛾(‖𝑥_𝑘−𝑥_𝑙‖²), and 𝛾^′_𝑘𝑙=𝛾^′(‖𝑥_𝑘−𝑥_𝑙‖²), (5) can be rewritten asd𝑥_𝑘=d𝑡^𝐿_𝑙=1𝛾_𝑘𝑙𝑎_𝑙,d𝑎_𝑘d𝑡=−2^𝐿_𝑙=1𝛾^′_𝑘𝑙(𝑎_𝑙⋅𝑎_𝑘)(𝑥_𝑘−𝑥_𝑙). (7) Once the initial position of the vertices, 𝑥(0)=(𝑥₁(0),…,𝑥_𝐿(0)), and the initial momentum, 𝑎(0)=(𝑎₁(0),…,𝑎_𝐿(0)), are provided, the evolution of the point set is uniquely determined.

3. Generative Model for Surface Observation

3.1. Random Triangulated Surfaces

If a triangulated template surface 𝑇 with vertices 𝑥^(𝑇) is given, and we solve, until time 𝑡=1, (5) initialized with 𝑥(0)=𝑥^(𝑇) and a random initial momentum 𝑎(0)=𝛼, the displaced vertices provide a random perturbation of the initial surface that will be denoted by 𝑇_𝛼. This is stated in the following definition.

Definition 1. Let 𝑇 be a triangulated surface with vertices 𝑥^(𝑇)=(𝑥₁^(𝑇),…,𝑥_𝐿^(𝑇)). Let 𝛼∈(ℝ³)^𝐿 be a collection of 𝐿 vectors in ℝ³. Let (𝑥(𝑡),𝑎(𝑡)) be the solution of (5) with initial condition 𝑥(0)=𝑥^(𝑇) and 𝑎(0)=𝛼. One defines 𝑇_𝛼 to be the triangulated surface with vertices 𝑥^(𝑇𝛼)=𝑥(1) and the same topology as 𝑇.

By letting 𝛼 be random, we build 𝑇_𝛼 as a random deformation of 𝑇. This will form the “ideal”, unobserved, surface, of which only a noisy version is observable (the noise process will be described in the next section).

Following [15], we will use a Bayesian formulation in which 𝑇 is itself represented as a random deformation 𝑇_0,𝜇∶=(𝑇₀)_𝜇, where 𝑇₀ is a fixed surface that we will call the hypertemplate, and 𝜇 is a prior initial momentum shooting from 𝑇₀ to 𝑇 (same notation as in Definition 1). One of the main interests of using a hypertemplate is to fix the topology of 𝑇 so that it belongs, by construction, to the same class of objects as 𝑇₀.

So, if 𝑁 surfaces are observed, we need to model the probability distribution of the prior momentum, 𝜇 (starting at 𝑇₀), which specifies 𝑇=𝑇_0,𝜇 and of 𝑁 deformation momenta 𝛼⁽¹⁾,…,𝛼^(𝑁) which specify the surfaces 𝑇_𝛼(1),…,𝑇_𝛼(𝑁). We now provide a statistical model for the joint probability distribution of 𝜇,𝛼⁽¹⁾,…,𝛼^(𝑁).

We first introduce some notation. Letting 𝐾 be the kernel introduced in the previous section to define the geodesic shooting equations, we let Γ_𝑇 be the 3𝐿 by 3𝐿 matrix formed with the 3 by 3 blocks 𝐾(𝑥_𝑘^(𝑇),𝑥_𝑙^(𝑇))𝐼𝑑_ℝ³. We define, for a triangulated surface 𝑇 with 𝐿 vertices 𝑥^(𝑇), and 𝛼∈ℝ^3𝐿, ‖𝛼‖²_𝑇=𝛼^∗Γ_𝑇𝛼=^𝐿_𝑘,𝑙=1𝐾(𝑥_𝑘^(𝑇),𝑥_𝑙^(𝑇))(𝛼_𝑘⋅𝛼_𝑙). (8)

We define the joint distribution of 𝜇,𝛼⁽¹⁾,…,𝛼^(𝑁) on ℝ^3𝐿×(ℝ^3𝐿)^𝑁 by𝑝(𝜇,𝛼⁽¹⁾,…,𝛼^(𝑁)1)=𝑍1exp(−2𝜆‖𝜇‖²_𝑇0−12^𝑁_𝑛=1‖𝛼^(𝑛)‖²_𝑇0,𝜇), (9) where 𝜆 is a fixed parameter regulating the weight on the hypertemplate.

There is a technical difficulty here, which is that one must make sure that this probability can be normalized (𝑍 exists), which requires that the exponential is integrable. That this is true is not straightforward, and we have not been able to find a proof that works with any choice of the kernel 𝐾. One way to deal with this is to introduce a constant 𝐴_𝜇 (which can be chosen arbitrarily large so that it does not interfere with the algorithms that will follow), and add to the model the constraint that ‖𝜇‖_𝑇0 is smaller than 𝐴_𝜇. Under such an assumption, one obtains (after integrating out the 𝛼's) 𝑍=(2𝜋)^3𝑁𝐿/2_{‖𝜇‖𝑇0≤𝐴𝜇}1exp(−2‖𝜇‖²_𝑇0−𝑁2logdetΓ_𝑇0,𝜇)𝑑𝜇. (10) This is finite, since, for any given 𝜇, the transformation 𝑥^(𝑇0)→𝑥^(𝑇0,𝜇) is the restriction of a diffeomorphism to the vertices of 𝑇₀ (as seen from (5)). This implies that Γ_𝑇0,𝜇 is nonsingular, and its determinant is bounded away from 0 when 𝜇 is restricted to a compact space.

In fact, the choice 𝐴_𝜇=∞ can be proved to be acceptable for a large class of kernels. Those are kernels for which the smallest eigenvalue of  Γ_𝑇 decreases at a speed which is at most polynomial in the minimal distance, ℎ_𝑇, between the vertices. A list of kernels satisfying this property can be found in [22]. For such kernels, we find that (𝐿 being fixed) logdetΓ_𝑇=𝑂(logℎ_𝑇). Just sketching the argument here, one can prove, using elementary properties of dynamical systems, that ℎ_𝑇=𝑂(exp(−𝐶‖𝜇‖_𝑇0,𝜇) for some constant 𝐶, so that the log determinant in (9) is linear in ‖𝜇‖_𝑇0,𝜇 and 𝑍 is well defined, even with 𝐴_𝑇0,𝜇=∞. For very smooth kernels, including the Gaussian, bounds on the smallest eigenvalue of Γ_𝑇 are much worse (with a decay which is exponential in (−1/ℎ²_𝑇)), and the previous argument does not work. Since the bounds in [22] hold uniformly with respect to the number of points, a polynomial bound may still be valid for a fixed 𝐿, although we were unable to discover it.

Notice that, conditionally to the template, the momenta 𝛼 are independent and follow a Gaussian distribution with inverse covariance matrix given by Γ_𝑇. An example of simulated random deformations obtained using such a model is provided in Figure 1.

(a)

(b)

Figure 1 

Random deformation of a template caudate surface (left). The six surfaces following the template are independent realization of the model described in (9).

3.2. Observation and Noise

The second part of our generative model is to describe the observation process, which takes an ideal surface 𝑇_𝛼 generated according to the model above, and returns the noisy observable.

Modeling such a noise process is a tricky issue. Obvious choices (like adding noise to the vertices of 𝑇_𝛼) do not work because one cannot assume that the observed surfaces are discretized consistently with the template. In this paper, we will work around this issue by assimilating the observation of a surface that of a singular Gaussian process.

For this, we consider that surfaces in ℝ³ are not observable directly, but only via their action on test functions, that we will call sensors. We define a sensor to be a smooth vector field 𝑤 over ℝ³ (typically with small support). Given an oriented surface 𝑆, define (𝑆,𝑤)=_𝑆𝑤(𝑠)⋅𝑁_𝑆(𝑠)𝑑𝜔_𝑆(𝑠), (11) where 𝑁_𝑆 is the normal to 𝑆. The real number, (𝑆,𝑤) is the measurement of 𝑆 through the sensor 𝑤.

Now, modeling noisy surfaces will result in assuming that, given any 𝑤, the measurement (𝑆,𝑤) is a random variable. We will assume that it is Gaussian, and more generally, that, given 𝑚 sensors, 𝑤₁,…,𝑤_𝑚, the random vector ((𝑆,𝑤_𝑗),𝑗=1,…,𝑚) is Gaussian.

𝑆, via its action on sensors, is therefore modeled as a Gaussian random field. Given an ideal surface 𝑇_𝛼, we will assume that its mean is given by𝐸((𝑆,𝑤))=(𝑇_𝛼,𝑤) (12) and the process is thereafter uniquely characterized by its covariance operator𝐺(𝑤,𝑤)=cov((𝑆,𝑤),(𝑆,𝑤)). (13) We will assume that this covariance is associated to a symmetric operator 𝐿_obs so that𝐺(𝑤,𝑤)=𝜎²_ℝ𝑑(𝐿_obs𝑤(𝑥)⋅𝑤(𝑥))𝑑𝑥𝑑𝑥. (14) (The apparently redundant parameter 𝜎², which could have been included in 𝐿_obs, appears here because it can be easily estimated by the algorithm, with the operator 𝐿_obs remaining constant.)

To finalize our model, it remains to describe how 𝑆 is discretized, that is, to make explicit a finite family of sensors through which 𝑆 is measured. Let (𝑧_𝑗,𝑗∈𝐽) form a regular grid of points in Ω. Let 𝛾_𝑠 be a radial basis function (a Gaussian, e.g.,) and define, for 𝑗∈𝐽 and 𝑑=1,2,3𝑤_𝑗,𝑑=𝛾_𝑠(𝑥−𝑧_𝑗)𝑒_𝑑 (15) where 𝑒_𝑑 is the 𝑑th vector of the standard basis of ℝ³ (this therefore specifies 3|𝐽| sensors). The resulting observed variables are𝑦_𝑗,𝑑=(𝑆,𝑤_𝑗,𝑑)=_𝑆𝛾_𝑠(𝑠−𝑧_𝑗)𝑁_𝑆^(𝑑)(𝑠)𝑑𝜔_𝑆(𝑠), (16) where 𝑁_𝑆^(𝑑)=𝑁_𝑆⋅𝑒_𝑑 is the 𝑑th coordinate of 𝑁_𝑆. These variables are, by assumption, jointly Gaussian, with means 𝑚_𝑗,𝑑=(𝑇_𝛼,𝑤_𝑗,𝑑) and covariance matrix𝑔_{𝑖,𝑑,𝑖′,𝑑′}=𝐺(𝑤_𝑖,𝑑,𝑤_{𝑖′,𝑑′})=𝜎²𝛿_{𝑑,𝑑′ℝ𝑑}𝐿_obs𝛾_𝑠(𝑥−𝑧_𝑖)𝛾_𝑠(𝑥−𝑧_𝑖′)𝑑𝑥𝑑𝑥. (17) Assuming that 𝐿_obs is translation invariant, the resulting expression is a function of 𝑧_𝑖−𝑧_𝑖′ that we will denote𝑔_{𝑖,𝑑,𝑖′,𝑑′}=𝜎²𝛿_𝑑𝑑′𝛾_obs(𝑧_𝑖−𝑧_𝑖′). (18)

Let 𝑅_obs=(𝑟^(obs)_𝑖𝑗) be the inverse matrix of the one with coefficients (𝛾_obs(𝑧_𝑗−𝑧^′_𝑗),𝑗,𝑗^′∈𝐽). The log likelihood of the process will include error terms taking the form ℰ_obs=1𝜎²_𝑗,𝑗′³_{∈𝐽𝑑=1}𝑟^(obs)_𝑗,𝑗′(𝑦_𝑗,𝑑−(𝑇_𝛼,𝑤_𝑗,𝑑))(𝑦_𝑗′,𝑑−(𝑇_𝛼,𝑤_𝑗′,𝑑)). (19) Replacing 𝑦_𝑗,𝑑 by its expression in (16), we have 𝜎²ℰ_obs=_𝑗,𝑗′³_{∈𝐽𝑑=1}𝑟^(obs)_{𝑗,𝑗′𝑆}𝛾_𝑠(𝑥−𝑧_𝑗)𝑁_𝑆^(𝑑)(𝑥)𝑑𝜔_𝑆×(𝑥)_𝑆𝛾_𝑠(𝑥−𝑧_𝑗′)𝑁_𝑆^(𝑑)(𝑥)𝑑𝜔_𝑆(𝑥)−2_{𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗,𝑗′³_𝑑=1𝑆𝛾_𝑠(𝑥−𝑧_𝑗)𝑁_𝑆^(𝑑)(𝑥)𝑑𝜔_𝑆×(𝑥)_𝑇𝛼𝛾_𝑠(𝑥−𝑧_𝑗′)𝑁_𝑇^(𝑑)_𝛼(𝑥)𝑑𝜔_𝑇𝛼(𝑥)×_{𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗𝑗′³_{𝑑=1𝑇𝛼}𝛾_𝑠(𝑥−𝑧_𝑗)𝑁_𝑇^(𝑑)_𝛼(𝑥)𝑑𝜔_𝑇𝛼(𝑥)×_𝑇𝛼𝛾_𝑠(𝑥−𝑧_𝑗′)𝑁_𝑇^(𝑑)_𝛼(𝑥)𝑑𝜔_𝑇𝛼(𝑥). (20) Let us analyze the first term. We have _{𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗,𝑗′³_𝑑=1𝑆𝛾_𝑠(𝑥−𝑧_𝑗)𝑁_𝑆^(𝑑)(𝑥)𝑑𝜔_𝑆×(𝑥)_𝑆𝛾_𝑠(𝑥−𝑧_𝑗′)𝑁_𝑆^(𝑑)(𝑥)𝑑𝜔_𝑆(𝑥)=_{𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗,𝑗′³_{𝑑=1𝑆×𝑆}𝛾_𝑠(𝑥−𝑧_𝑗)𝑁_𝑆^(𝑑)(𝑥)𝛾_𝑠×(𝑦−𝑧_𝑗′)𝑁_𝑆^(𝑑)(𝑦)𝑑𝜔_𝑆(𝑥)𝑑𝜔_𝑆(𝑦)=_{𝑆×𝑆𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗,𝑗′𝛾_𝑠(𝑥−𝑧_𝑗)𝛾_𝑠(𝑦−𝑧_𝑗′)×(𝑁_𝑆(𝑥)⋅𝑁_𝑆(𝑦))𝑑𝜔_𝑆(𝑥)𝑑𝜔_𝑆(𝑦)=_𝑆×𝑆𝐾_obs(𝑥,𝑦)(𝑁_𝑆(𝑥)⋅𝑁_𝑆(𝑦))𝑑𝜔_𝑆(𝑥)𝑑𝜔_𝑆(𝑦), (21) with the notation𝐾_obs(𝑥,𝑦)=_{𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗,𝑗′𝛾_𝑠(𝑥−𝑧_𝑗)𝛾_s(𝑦−𝑧_𝑗′). (22)

Treating the other terms similarly, we can rewrite the error term in the form ℰ_obs=1𝜎²_𝑆×𝑆𝐾_obs(𝑥,𝑦)(𝑁_𝑆(𝑥)⋅𝑁_𝑆(𝑦))𝑑𝜔_𝑆(𝑥)𝑑𝜔_𝑆−2(𝑦)𝜎²_{𝑆×𝑇𝛼}𝐾_obs(𝑥,𝑦)×(𝑁_𝑆(𝑥)⋅𝑁_𝑇𝛼(𝑦))𝑑𝜔_𝑆(𝑥)𝑑𝜔_𝑇𝛼+1(𝑦)𝜎²_{𝑇𝛼×𝑇𝛼}𝐾_obs(𝑥,𝑦)×(𝑁_𝑇𝛼(𝑥)⋅𝑁_𝑇𝛼(𝑦))𝑑𝜔_𝑇𝛼(𝑥)𝑑𝜔_𝑇𝛼(𝑦) (23) and we will abbreviate this (introducing a notation for the right-hand side) as ℰ_obs=(1/𝜎²)‖𝑆−𝑇_𝛼‖²_obs. Thus, we can write 𝑝(𝑦∣𝑇_𝛼1,𝜎)=cstexp(−2𝜎²‖𝑆−𝑇_𝛼‖²_obs). (24)

We have the following important proposition.

Proposition 1. Assume that 𝛾_𝑠 is taken equal to the Green function of 𝐿_𝑜𝑏𝑠. Then, when the grid 𝐽 becomes finer and 𝐾_𝑜𝑏𝑠 is given by (22), one has 𝐾_obs(𝑥,𝑦)→𝛾_𝑜𝑏𝑠(𝑥−𝑦). (25)

See the appendix for a proof of this proposition (which requires some background on the theory of Hilbert spaces with a reproducing kernel) and possible extensions. For practical purposes, we will use 𝛾_obs instead of 𝐾_obs in ‖⋅‖_obs, therefore assuming that the sensors are chosen according to the proposition. It is interesting to notice that the resulting norm in this case is precisely the norm that has been introduced in [16] to compare surfaces, when they are considered as elements of a reproducing kernel Hilbert space of currents. We refer to [16] for details on the mathematical construction.

In the rest of the paper, we develop a parametric procedure to estimate the template 𝑇 from the observation of i.i.d. surfaces 𝑆⁽¹⁾,𝑆⁽²⁾,…,𝑆^(𝑁) generated as described above. This includes in particular 𝑁 hidden deformation momenta 𝛼⁽¹⁾,𝛼⁽²⁾,…,𝛼^(𝑁), such that the complete distribution of observed and unobserved variables is𝑝(𝜇,𝛼⁽¹⁾,𝑆⁽¹⁾,…,𝛼^(𝑁),𝑆^(𝑁))𝜆=cstexp(−2‖𝜇‖²_𝑇0−12^𝑁_𝑛=1(‖𝛼^(𝑛)‖²_𝑇+1𝜎²‖𝑆^(𝑛)−𝑇_𝛼(𝑛)‖²_obs)), (26) where we have written for short 𝑇=𝑇_0,𝜇.

We now discuss the estimation of the parameter 𝜇, and of the associated template 𝑇_0,𝜇, which is the main purpose of this paper. This will be implemented with a mode approximation of the EM algorithm, as described in the next section.

4. Algorithms for Surface Template Estimation

4.1. Mode Approximation EM

Let Θ=(𝛼⁽¹⁾,𝛼⁽²⁾,…,𝛼^(𝑁)) be the hidden part of the process, representing the collection of initial momenta, and let 𝐒=(𝑆⁽¹⁾,𝑆⁽²⁾,…,𝑆^(𝑁)) be the collection of observed surfaces. The complete distribution for the process, including the prior is given by (26).

The EM algorithm is an iterative method that updates a current estimation of 𝜇 using the following two 𝐸- and 𝑀-steps.

𝐸-step: determine the conditional expectation 𝜇↦𝐸_𝜇{log𝜋(𝜇,Θ,𝐒)∣𝐒)}.

𝑀-step: maximize this expression with respect to 𝜇 and replace the current estimation of 𝜇 by the obtained maximizer.

The conditional expectation can be expanded as𝐸_𝜇𝜆{log𝜋(𝜇,Θ,𝐒)∣𝐒)}=−2‖𝜇‖²_𝑇0−12^𝑁_𝑛=1𝐸_𝜇(‖𝛼^(𝑛)‖²_𝑇+1𝜎²‖𝑆^(𝑛)−𝑇_𝛼(𝑛)‖²_obs∣𝑆^(𝑛)). (27)

Considering the highly nonlinear relation between 𝛼^(𝑛) and the deformed surface 𝑇_𝛼(𝑛), an explicit computation in (27) is impossible. We will therefore rely on the classic mode approximation in the EM which replaces the conditional distribution by a Dirac measure taken at the conditional mode, yielding𝐸_𝜇≈𝜆{log𝜋(𝜇,Θ,𝐒)∣𝐒)}2‖𝜇‖²_𝑇0+12^𝑁_𝑛=1‖𝛼^(𝑛)‖²_𝑇+1𝜎²‖𝑆^(𝑛)−𝑇𝛼^(𝑛)‖²_obs}, (28) where 𝛼^(𝑛)=argmin_𝛼{‖𝛼‖²_𝑇+1𝜎²‖𝑆^(𝑛)−𝑇_𝛼‖²_obs}. (29)

This results in a maximum a posteriori procedure that maximizes alternatively in 𝜇 and in the 𝛼^(𝑛)'s. Like in [15], we will refer to it as a mode approximation to the EM (MAEM) rather than a MAP algorithm, in order to strengthen the fact that it is an approximation of the maximum a posteriori procedure relying on the likelihood of the observed data only. As illustrated in [14], the MAEM can be biased (leading to inexact estimation of the template, even with a very large number of samples), especially when the noise is important, but it is obviously more feasible than the exact EM. Notice that the MAEM method is a special form of EM algorithm, and as such optimizes a lower bound of the log likelihood of the observed data.

We summarize the two steps of the 𝑖th iteration of the MAEM in our case. Suppose 𝜇 and the 𝛼^(𝑛)'s are the current variables to be updated. Then the next iteration is as follows.

MAE step: with 𝜇 (and therefore 𝑇) fixed, find, for 𝑛=1,…,𝑁, 𝛼^(𝑛) to minimize‖𝛼‖²_𝑇+1𝜎²‖𝑆^(𝑛)−𝑇_𝛼‖²_obs (30) and replace 𝛼^(𝑛) by 𝛼^(𝑛).

𝑀 step: with 𝛼^(𝑛) fixed, update 𝜇 with the minimizer of𝜆‖𝜇‖²_𝑇0+^𝑁_𝑛=1(‖𝛼^(𝑛)‖²_𝑇0,𝜇+1𝜎²‖𝑆^(𝑛)−(𝑇₀𝜇)_𝛼(𝑛)‖²_obs). (31)

We now discuss how each of these steps can be implemented.

4.2. MAE Step

Our goal in this section is to optimize (30). We work with fixed 𝑛 and drop it from the notation to simplify the expressions. The objective function is𝐸(𝛼)=‖𝛼‖²_𝑇+1𝜎²‖𝑆−𝑇_𝛼‖²_obs. (32) This problem is equivalent to the surface matching algorithm considered in [16], with a slightly different formulation since [16] optimize an energy with respect to a time-dependent momentum instead of just the initial momentum (i.e., they solved simultaneously the geodesic estimation and the matching problems). These two formulations are equivalent when using continuous time (they produce the same minima), but they yield different results when discretized. In our setting, formulating the problem as in (32) is natural, and focuses on the modeled random variable, 𝛼.

We need to compute the variation of 𝐸 with respect to 𝛼. This computation will be useful for the 𝑀-step also. We first discuss the discretization of the error term, which follows [16]. Let 𝑆 be a triangulated surfaces, with vertices 𝑥^(𝑆)=(𝑥₁^(𝑆),…,𝑥_𝐿^(𝑆)) and faces 𝐹^(𝑆)=(𝑓₁^(𝑆),…,𝑓_𝑀^(𝑆)). Each face is represented by an ordered triple of vertices: 𝑓=(𝑥^𝑓₁,𝑥^𝑓₂,𝑥^𝑓₃), and we define the face centers and area-weighted normals by 𝑐_𝑓=13(𝑥^𝑓₁+𝑥^𝑓₂+𝑥^𝑓₃𝑁),_𝑓=12(𝑥^𝑓₂−𝑥^𝑓₁)×(𝑥^𝑓₃−𝑥^𝑓₁). (33) Then a discrete approximation of ‖𝑆−𝑆^′‖²_obs is 𝑈(𝑥^𝑆′)=_{𝑓,𝑓′∈𝐹(𝑆)}𝛾_obs(𝑐_𝑓−𝑐_𝑓′)(𝑁_𝑓⋅𝑁_𝑓′)+_{𝑓,𝑓′∈𝐹′)(𝑆}𝛾_obs(𝑐_𝑓−𝑐_𝑓′)(𝑁_𝑓⋅𝑁_𝑓′)−2_{𝑓∈𝐹(𝑆),𝑓′∈𝐹′)(𝑆}𝛾_obs(𝑐_𝑓−𝑐_𝑓′)(𝑁_𝑓⋅𝑁_𝑓′), (34) where 𝑆 is considered as fixed and 𝑈 is therefore considered as a function of the vertices, 𝑥^(𝑆′), of 𝑆^′.

With this notation, we can write 𝐸(𝛼)=‖𝛼‖²_𝑇+1𝜎²𝑈(𝑥^(𝑇𝛼)). (35) We want to compute the gradient of 𝐸 and for this, apply the chain rule to the transformations 𝛼→𝑇_𝛼→𝑈(𝑥^(𝑇𝛼)), yielding𝜕_𝛼𝑈=(𝑑𝑥^(𝑇𝛼))𝑑𝛼^∗𝜕_𝑥𝑈, (36) where 𝐴^∗ is the transpose matrix of 𝐴.

The gradient of 𝑈 with respect to 𝑥^(𝑆′) has been computed in [16], and is given as follows. We denote by 𝐹_𝑙^(𝑆) the set of faces (triangles) that contain a vertex 𝑥_𝑙 in a triangulated surface 𝑆. For 𝑓∈𝐹_𝑙^(𝑆), we let 𝑒_𝑙(𝑓) denote the edge opposed to 𝑥_𝑙 in 𝑓, positively oriented in 𝑓. With 𝑆^′=𝑇_𝛼, we have𝜕𝑈𝜕𝑥^(𝑆′)_𝑙(𝑥^(𝑆′))=_{𝑓∈𝐹′)𝑙(𝑆𝑔∈𝐹(𝑆)∪𝐹′)(𝑆}𝑎_𝑔(23∇𝛾_obs(𝑐_𝑓−𝑐_𝑔)(𝑁_𝑓⋅𝑁_𝑔)+𝛾_obs(𝑐_𝑓−𝑐_𝑔)𝑒_𝑙(𝑓)×𝑁_𝑔), (37) where 𝑎_𝑔=1 if 𝑔∈𝐹^(𝑆′) and 𝑎(𝑔)=−1 if 𝑔∈𝐹^(𝑆).

Now let us derive the variation of 𝑥^(𝑇𝛼) with respect to the initial momentum, 𝛼. We know that 𝑥^(𝑇𝛼)=𝑥(1), where 𝑥 and 𝑎 evolve according to the system (7) d𝑥_𝑘=d𝑡^𝐿_𝑙=1𝛾_𝑘𝑙𝑎_𝑙,d𝑎_𝑘d𝑡=−2^𝐿_𝑙=1𝛾^′_𝑘𝑙(𝑎_𝑙⋅𝑎_𝑘)(𝑥_𝑘−𝑥_𝑙) (38) with 𝑥(0)=𝑥^(𝑇) and 𝑎(0)=𝛼 (and 𝛾_𝑘𝑙,𝛾^′_𝑘𝑙 are short for 𝛾(‖𝑥_𝑘−𝑥_𝑙‖²) and 𝛾^′(‖𝑥_𝑘−𝑥_𝑙‖²)). Now an infinitesimal variation 𝛼→𝛼+𝛿𝛼 in the initial condition induces infinitesimal variations 𝑎+𝛿𝑎 and 𝑥+𝛿𝑥 over time, and the pair (𝛿𝑥,𝛿𝑎) obeys the following differential system, that can be obtained from a formal differentiation of (7):d𝛿𝑥_𝑘(𝑡)=d𝑡^𝐿_𝑙=1𝛾_𝑘𝑙𝛿𝑎_𝑙+2^𝐿_𝑙=1𝛾^′_𝑘𝑙𝑎_𝑙(𝑥_𝑘−𝑥_𝑙)⋅(𝛿𝑥_𝑘−𝛿𝑥_𝑙),d𝛿𝑎_𝑘(𝑡)d𝑡=−2^𝐿_𝑙=1𝛾^′_𝑘𝑙(𝑎_𝑙⋅𝛿𝑎_𝑘+𝛿𝑎_𝑙⋅𝑎_𝑘)(𝑥_𝑘−𝑥_𝑙)−2^𝐿_𝑙=1𝛾^′_𝑘𝑙(𝑎_𝑙⋅𝑎_𝑘)(𝛿𝑥_𝑘−𝛿𝑥_𝑙)−4^𝐿_𝑙=1𝛾^′′_𝑘𝑙(𝑎_𝑙⋅𝑎_𝑘)(𝑥_𝑘−𝑥_𝑙)((𝑥_𝑘−𝑥_𝑙)⋅(𝛿𝑥_𝑘−𝛿𝑥_𝑙)) (39)(40) with 𝛾^′′_𝑘𝑙=𝛾^′′(‖𝑥_𝑘−𝑥_𝑙‖²).

One can rewrite it in the matrix form:dd𝑡𝛿𝑥𝛿𝛼=𝐽(𝑡)𝛿𝑥𝛿𝛼, (41) where 𝐽(𝑡)=(^{𝐽𝑥𝑥𝐽𝑥𝑎}_{𝐽𝑎𝑥𝐽𝑎𝑎}) with𝐽_𝑥𝑥(𝑘,𝑙)=(2^𝐿_𝑞=1𝛾^′_𝑘𝑞𝑎_𝑞(𝑥_𝑙−𝑥_𝑞)^∗)𝛿_𝑘𝑙−2𝛾^′_𝑘𝑙𝑎_𝑙(𝑥_𝑘−𝑥_𝑙)^∗𝐽_𝑥𝑎(𝑘,𝑙)=𝛾_𝑘𝑙𝐼𝑑_ℝ³.𝐽_𝑎𝑥(𝑘,𝑙)=−(2^𝐿_𝑞=1(𝑎_𝑞⋅𝑎_𝑙)×(𝛾^′_𝑘𝑞𝐼𝑑_ℝ³+2𝛾^′′_𝑘𝑞(𝑥_𝑙−𝑥_𝑞)(𝑥_𝑙−𝑥_𝑞)^∗)𝛿_𝑘𝑙)+2(𝑎_𝑙⋅𝑎_𝑘)(𝛾^′_𝑘𝑙𝐼𝑑_ℝ³+2𝛾^′′_𝑘𝑙(𝑥_𝑘−𝑥_𝑙)(𝑥_𝑘−𝑥_𝑙)^∗𝐽)._𝑎𝑎(𝑘,𝑙)=−(2(^𝐿_𝑞=1𝛾^′_𝑘𝑞(𝑥_𝑙−𝑥_𝑞)𝑎^∗_𝑞))𝛿_𝑘𝑙−2𝛾^′_𝑘𝑙(𝑥_𝑘−𝑥_𝑙)𝑎^∗_𝑘. (42) Solving this system with initial condition 𝛿𝑥(0)=0 and 𝛿𝑎(0)=𝛿𝛼 provides what we have denoted (𝑑𝑥^(𝑇𝛼)d𝛼)𝛿𝛼. (43) One does not need to compute all the coefficients of the matrix 𝑑𝑥^(𝑇𝛼)/𝑑𝛼 using this equation in order to apply the transpose in (36). This is fortunate because this would constitute a computationally demanding effort given that this matrix is 3𝐿 by 3𝐿 with 𝐿 large. The right hand side of (36) can be in fact computed directly using a single dynamical system, given byd𝜂d𝑡𝜂^𝑥_𝛼=−𝐽(𝑡)^∗𝜂_𝑥𝜂_𝛼, (44) where 𝐽(𝑡) is defined in (39). If (44) is solved from time 𝑡=1 to time 𝑡=0 with 𝜂_𝑥(1)=𝜕_𝑥𝑈 and 𝜂_𝛼(1)=0, then 𝜕_𝛼𝑈=𝜂_𝛼(0). (45) This is a simple consequence of the theory of linear differential systems (a proof is provided in the Appendix for completeness). Note that the matrix 𝐽(𝑡) depends on the solution of (7) computed with initial conditions 𝑥(0)=𝑥^(𝑇) and 𝑎(0)=𝛼. To emphasize this dependency, we will denote it 𝐽(𝑡)=𝐽^(𝑇,𝛼)(𝑡) in the following.

Given this, we see that a variation 𝛼→𝛼+𝛿𝛼 induces a first-order variation 𝛿𝐸 of the energy given by𝛿𝐸=2⟨𝛿𝛼,𝛼⟩_𝑇+1𝜎²𝛿𝛼⋅𝜂_𝛼(0), (46) where the 𝑇-dot product is ⟨𝛿𝛼,𝛼⟩_𝑇=𝛿𝛼⋅(Γ_𝑇𝛼) and Γ_𝑇 is the matrix formed with 3 by 3 blocs 𝛾_𝑘𝑙𝐼𝑑_ℝ³.

We choose to operate the gradient descent with respect to this dot product and therefore choose a variation proportional to 𝛿𝛼=−(2𝛼+Γ_𝑇⁻¹𝜂_𝛼(0)). So, the algorithm to compute an optimal 𝛼 is the following.

Algorithm 1 (MAE Step for Surface Template Estimation). (1)Compute the variation 𝜕_{𝑥𝛼)(𝑇}𝑈 using (37). (2)Solve backward in time (44) initialized with 𝜂_𝑥(1)=𝜕_{𝑥𝛼)(𝑇}𝑈 and 𝜂_𝛼=0.(3)Replace 𝛼 by 𝛼−𝜖(2𝛼+(1/𝜎²)Γ_𝑇⁻¹𝜂_𝛼(0)) (using a line-search to optimize 𝜖).

This algorithm has to be applied 𝑁 times (for all 𝛼_𝑘, 𝑘=1,…,𝑁) in the MAE step.

Remark 1. The matrix Γ_𝑇 being typically very badly conditioned, we numerically compute Γ_𝑇⁻¹𝜂_𝛼 after adding a small positive number to the diagonal of Γ_𝑇. The inversion itself is computed using conjugate gradient.

4.3. M Step

There are many similarities between the 𝑀-step and the 𝐸-step variational problems, so that we will be able to only sketch the detail of the computation here. We need to minimize 𝐸(𝜇)=‖𝜇‖²_𝑇0+1𝜆^𝑁_𝑛=1𝑈^(𝑛)(𝑥^(𝑇0,𝜇)) (47) with𝑈^(𝑛)(𝑥^(𝑇))=‖𝛼^(𝑛)‖²_𝑇+1𝜎²‖𝑆−𝑇_𝛼(𝑛)‖²_obs. (48)

Let us consider the variation of each term in the sum (fixing 𝑛, that we temporarily drop from the notation). Since ‖𝛼‖²_𝑇=^𝐿_𝑘,𝑙=1𝛾(‖𝑥_𝑘^(𝑇)−𝑥_𝑙^(𝑇)‖²)(𝛼_𝑘⋅𝛼_𝑙), (49) we can write𝜕‖𝛼‖²_𝑇𝜕𝑥_𝑘^(𝑇)=2^𝐿_𝑙=1𝛾^′(‖𝑥_𝑘^(𝑇)−𝑥_𝑙^(𝑇)‖²)(𝛼_𝑘⋅𝛼_𝑙)(𝑥_𝑘^(𝑇)−𝑥_𝑙^(𝑇)). (50)

The function 𝑈 being defined as before, we see that the derivative of the second term is given by applying the chain rule again, this time in the form (𝑑𝑥^(𝑇𝛼)𝑑𝑥^(𝑇))^∗𝜕_{𝑥𝛼)(𝑇}𝑈. (51)

Like in the previous section, the transpose of the differential applied to the gradient of 𝑈 can be computed by solving a dynamical system backward in time. In fact, it is the same system as with the variation in 𝛼, namely, d𝜂d𝑡𝜂^𝑥_𝛼=−𝐽^(𝑇,𝛼)(𝑡)^∗𝜂_𝑥𝜂_𝛼 (52) still initialized with 𝜂_𝑥(1)=𝜕_{𝑥𝛼)(𝑇}𝑈 and 𝜂_𝛼(1)=0, but the relevant result now is 𝜂_𝑥(0). The gradient of 𝑈 is then𝜕_𝑥(𝑇)𝑈=𝜕_𝑥(𝑇)‖𝛼‖²_𝑇+1𝜎²𝜂_𝑥(0). (53)

Once this is computed, the next step is to compute (reintroducing 𝑛 in the notation) (𝑑𝑥^(𝑇0𝜇))𝑑𝜇^∗(^𝑁_𝑛=1𝜕_𝑥(𝑇)𝑈^(𝑛)). (54) This follows a similar procedure, using (44), with 𝑇₀ instead of 𝑇 and 𝜇 instead of 𝛼. This requires solvingd𝜂d𝑡𝜂^𝑥_𝜇=−𝐽^(𝑇0,𝜇)(𝑡)^∗𝜂_𝑥𝜂_𝜇, (55) initialized with 𝜂_𝑥(1)=𝜕_𝑥(𝑇)𝑈 and 𝜂_𝜇(1)=0. The variation of 𝐸 associated to an infinitesimal variation of 𝜇 is then 𝛿𝐸=⟨𝛿𝜇,2𝜇⟩_𝑇0+1𝜆(𝛿𝜇⋅𝜂_𝜇1(0))=⟨𝛿𝜇,2𝜇+𝜆Γ_𝑇⁻¹₀𝜂_𝜇(0)⟩. (56)

We summarize the M step in the following algorithm.

Algorithm 2 (𝑀-Step Algorithm for Surface Template Estimation). (1)For 𝑛=1,…,𝑁: (1.1)Compute 𝜕_𝑥𝛼^)(𝑇_(𝑛)𝑈 using (37). (1.2)Solve system (44) backward in time with initial condition 𝜂_𝑥^(𝑛)(1)=𝜕_𝑥𝛼^)(𝑇_(𝑛)𝑈 and 𝜂_𝛼^(𝑛)(1)=0. (1.3)Compute 𝜕_𝑥(𝑇)𝑈^(𝑛)=𝜕_𝑥(𝑇)‖𝛼^(𝑛)‖²_𝑇+1𝜎²𝜂_𝑥^(𝑛)(0) (57) using (50).(2)Solve system (55) backward in time with 𝜂_𝑥(1)=^𝑁_𝑛=1𝜕_𝑥(𝑇)𝑈^(𝑛) (58) and 𝜂_𝜇(1)=0.(3)Replace 𝜇 by 𝜇−𝜖(𝜇+Γ_𝑇⁻¹₀𝜂_𝜇(0)/𝜆), 𝜖 being optimized with a line search.

4.4. Surface Template Estimation Algorithm

We finally summarize the surface template estimation algorithm:

Algorithm 3 (Surface Template Estimation). Having the hypertemplate 𝑇₀ and observed surfaces 𝑆⁽¹⁾,…,𝑆^(𝑁), the goal is to estimate the template 𝑇. Let 𝑇,𝜇,𝛼⁽¹⁾,…,𝛼^(𝑁) denote the current estimation with initial guess 𝑇=𝑇₀, 𝜇=0, 𝛼^(𝑛)=0. Then, in the next iteration, update with the following steps: (i)With 𝑇 fixed, apply Algorithm 1 to update each 𝛼^(𝑛), 𝑛=1,…,𝑁. (ii)With 𝛼^(𝑛)'s fixed, apply Algorithm 2 to obtain a new value for 𝜇. (iii)Solve (7) initialized with the hypertemplate 𝑇₀ and the newly obtained 𝜇 to update the estimated template, 𝑇.

5. Result and Discussion

We applied the algorithm to surface data of human brain's caudate, thalamus, and hippocampus. All data are courtesy of Center for Imaging Science at Johns Hopkins University. Each surface has around 5–10 thousand triangle cells. We randomly chose one as the hypertemplate and the others as observed surfaces. In these experiments, we set 𝜆=1.0 and 𝜎²=1.0.

Figures 2 and 3 are the template estimation result for caudate and thalamus, respectively.

(a) 𝑆1

(b) 𝑆2

(c) 𝑆3

(d) 𝑆4

(e) 𝑆5

(f) 𝑆6

(g) 𝑆7

(h) 𝑆8

(i) 𝑆9

(j) Hypertemplate 𝑇0

(k) Estimated template 𝑇

Figure 2 

Estimating the surface template from 9 caudate data. (a)–(i) observed surfaces. (j) is the hypertemplate. (k) is the result.

(a) 𝑆1

(b) 𝑆2

(c) 𝑆3

(d) 𝑆4

(e) 𝑆5

(f) 𝑆6

(g) 𝑆7

(h) 𝑆8

(i) 𝑆9

(j) Hypertemplate 𝑇0

(k) Estimated template 𝑇

Figure 3 

Estimating the surface template from 9 thalamus data. (a)–(i) observed surfaces. (j) is the hypertemplate. (k) is the result.

We also applied the algorithm to 101 hippocampus surface data in the BIRN Project (Biomedical Informatics Research Network). In Figure 4, Panels (a)–(h) are 8 examples of the 101 observations. Panel (i) is the hypertemplate and Panel (j) is the estimated template.

(a) 𝑆1

(b) 𝑆2

(c) 𝑆3

(d) 𝑆4

(e) 𝑆5

(f) 𝑆6

(g) 𝑆7

(h) 𝑆8

(i) Hypertemplate 𝑇0

(j) Estimated template 𝑇

Figure 4 

Estimating the surface template from 101 data. (a)–(h) are 8 examples out of 101 observed surfaces. (i) is the hypertemplate. (j) is the result.

The result is visually satisfying in the sense that the estimated template is found to agree with a qualitative representation of the population. For example, in the caudate experiment, the estimated template has an obviously narrower upper tip than the hypertemplate. This captures the population characteristic since most observed data have narrower upper tips. Notice that the obtained template does not look smoother than the rest of the population, as would typically yield template estimation methods that average over fixed parametrizations.

Figure 5 shows how the energy in (31) changes with the iteration in the caudate experiment. This confirms the effectiveness and convergence of the algorithm. One can see the energy drops quickly in the first twenty iterations, then gradually slows down. After the 35th iteration, the energy changes little and the estimated template remains stable.

Figure 5 

Energy change with the iteration.

In our model, the hypertemplate can be provided by an atlas obtained from other studies, although we here simply choose one of the surfaces in the population. Actually, as Figure 6 shows, different choices for the hypertemplate yield very similar results.

(a) Hypertemplate 1

(b) Result 1

(c) Hypertemplate 2

(d) Result 2

(e) Hypertemplate 3

(f) Result 3

Figure 6 

For the same observed population, we choose different surfaces as hypertemplate. The results only have minor differences.

6. Conclusion

In this paper, we have presented a Bayesian approach for surface template estimation. We have built, for this purpose, a generative statistical model for surfaces: the construction first applies a random initial momentum to the template surface, then assumes an observation process using test functions and noise. The template is assumed to be generated as a random deformation of a hypertemplate, completing the Bayesian model. We used a mode approximation EM scheme to estimate the surface template, and introduced for this purpose a novel surface matching algorithm optimizing with respect to the initial momentum. The procedure has been tested with caudate, thalamus, and hippocampus surface data, showing its effectiveness and convergence, and also experimentally proved to be robust to variations in the choice of the hypertemplate.

Appendices

A. Proof of Proposition 1

Let us assume that 𝐿_obs𝛾_𝑠=𝛿₀, that is, 𝛾_𝑠 is the Green Kernel of the operator 𝐿_obs. This assumption implies, in particular, that 𝛾_𝑠=𝛾_obs. Given a smooth function 𝑓, denote 𝑓_𝐽(𝑥)=_{𝑗,𝑗′∈𝐽}𝑟^(obs)_𝑗,𝑗′𝛾_obs(𝑥−𝑧_𝑗)𝑓(𝑧_𝑗′). (A.1) Define the vector space 𝑉_𝐽={𝑔∈𝑉∶𝑔(𝑥)=_𝑗∈𝐽𝛾_obs(𝑥−𝑧_𝑗)𝜌_𝑗,𝜌_𝑗∈ℝ,𝑗∈𝐽}. (A.2) Introduce the RKHS 𝑉 with scalar product given by⟨𝑤,𝑤^′⟩_𝑉=_ℝ𝑑(𝐿_obs𝑤)𝑤^′𝑑𝑥. (A.3) Then 𝑓_𝐽 can be identified to the orthogonal projection of 𝑓 on 𝑉_𝑗, because 𝑟^(obs) is the inverse of the Gram matrix (matrix of dot products) of (𝛾_obs(⋅−𝑧_𝑗),𝑗∈𝐽), which are the generators of 𝑉_𝐽, and𝑓(𝑧_𝑗)=⟨𝑓,𝛾_obs(⋅,𝑧_𝑗)⟩_𝑉 (A.4) by assumption.

Now, let 𝐽_𝑚 be a sequence of progressively finer grids with resolution 𝜖_𝑚 tending to 0 when 𝑚 tends to infinity. We prove that ‖𝑓_𝐽𝑚−𝑓‖_𝑉→0, for which it suffices to prove that ⋃𝑉_𝐽𝑚 is dense in 𝑉. This is equivalent to the fact that no nonzero 𝑔 in 𝑉 can be orthogonal to all the 𝑉_𝐽𝑚's. But 𝑔 orthogonal to 𝑉_𝐽𝑚 is only possible when, for all 𝑗∈𝐽_𝑚, 𝑔(𝑧_𝑗)=⟨𝑔,𝛾_𝑤(⋅−𝑧_𝑗)⟩_𝑉=0. (A.5) Since the 𝑧^′_𝑗𝑠 form an arbitrarily fine grid in 𝑉 and functions in 𝑉 are continuous, this implies 𝑔=0.

So, 𝑓_𝐽→𝑓 in 𝑉, which implies, for example, pointwise convergence as soon as 𝑉 is embedded in the set of continuous functions (that is, if 𝛾_obs is continuous). This directly implies Proposition 1 in the case 𝛾_𝑠=𝛾_obs, by taking 𝑓(𝑥)=𝛾_obs(𝑥−𝑦) for a given 𝑦.

Extensions of this result is when 𝛾_𝑠 is obtained by applying some linear operator, say 𝐴, to 𝛾_obs. One then has𝐾_obs(⋅,𝑥)=𝐴(_{𝑗,𝑗′∈𝐽}𝑟^(𝑤)_𝑗,𝑗′𝛾_obs(⋅−𝑧_𝑗)𝛾_𝑠(𝑥−𝑧_𝑗′)) (A.6) and passing to the limit in the sum (for which one needs 𝛾_𝑠∈𝑉 and 𝐴 continuous on 𝑉), one gets 𝐾_obs(𝑥,𝑥)→𝐴²𝛾_obs.

B. Transposing Linear Differential Equations

We here justify the procedure described in Section 4, and prove the following fact: consider the solution, 𝑧(𝑡)∈ℝ^𝑝, of the ordinary differential equation 𝜕_𝑡𝑧=𝐽(𝑡)𝑧, (B.1) where  𝐽 is a time-dependent operator (a 𝑝 by 𝑝 matrix). Then, for any 𝑢∈ℝ^𝑝, we have 𝑢⋅𝑧(1)=𝜂(0)⋅𝑧(0), (B.2) where 𝜂 is the solution of the differential equation 𝜕_𝑠𝜂=−𝐽(𝑠)^∗𝜂 (B.3) with 𝜂(1)=𝑢.

Let us prove this result. First remark that since 𝑧 is the solution of a linear equation, it depends linearly on the initial condition, 𝑧(0). More precisely, let 𝑀(𝑠,𝑡) be the 𝑝 by 𝑝 matrix such that 𝜕_𝑡𝑀(𝑠,𝑡)=𝐽(𝑡)𝑀(𝑠,𝑡) (B.4) and 𝑀(𝑠,𝑠)=𝐼𝑑_ℝ^𝑝. Then 𝑧(𝑡)=𝑀(0,𝑡)𝑧(0), and, obviously, 𝜂(0)=𝑀(0,1)^∗𝑢. Using the identity 𝑀(𝑡,𝑠)𝑀(𝑠,𝑡)=𝐼𝑑_ℝ^𝑝, we obtain 0=𝑀(𝑡,𝑠)(𝜕_𝑠𝑀(𝑠,𝑡))+(𝜕_𝑠𝑀(𝑡,𝑠))𝑀(𝑠,𝑡)=(𝜕_𝑠𝑀(𝑠,𝑡))𝑀(𝑡,𝑠)+𝐽(𝑠) (B.5) so that 𝜕_𝑠𝑀(𝑠,𝑡)=−𝑀(𝑠,𝑡)𝐽(𝑠). Computing the transposed equation yields and taking 𝑡=1𝜕_𝑠𝑀(𝑠,1)^∗=−𝐽(𝑠)^∗𝑀(𝑠,1)^∗. (B.6) Thus, if 𝜂(𝑠)=𝑀(𝑠,1)^∗𝑢, we have 𝜂(1)=𝑢 and 𝜕_𝑠𝜂=−𝐽(𝑠)^∗𝜂 as announced.

Acknowledgments

This paper is partially supported by R01-MH056584, R01-EB000975, R01-MH60883, P50-MH071616, P01-AG003991, P01-AG0262761, P41-RR015241, R24-RR021382, and NSF DMS-0456253.