Finite dimensional landmark configurations are essential in shape analysis and computational anatomy, both for marking and following anatomically important areas in e.g. changing brain anatomies and discretely represented curve outlines, and in being among the simplest non-linear shape spaces. This simplicity, in particular the finite dimensionality, makes landmarks useful for theoretical investigations and for deriving algorithms that can subsequently be generalized to infinite dimensional shape spaces.
While probability distributions in Euclidean space can often be specified conveniently from their density function, e.g. the normal distribution with the density, the non-linear nature of shape spaces often rules out closed form functions. Indeed, a density defined in coordinates will be dependent on the chosen coordinate chart and thus not geometrically intrinsic, and normalization factors can be inherently hard to compute. A different approach defines probability distributions as transition distributions of stochastic processes. Because stochastic differential equations (SDEs) can be specified locally from their infinitesimal variations, it is natural to define them in geometric spaces. Belonging to this category, the present paper aligns with a range of recent research activities on nonlinear SDEs in shape analysis and geometric mechanics [26, 25, 14, 3, 4].
We consider here observations distributed according to the transition distribution of a Brownian motion, which is arguably one of the most direct generalizations of the Gaussian distribution to nonlinear geometries. For the Brownian motion, each infinitesimal step defining the SDE can be considered normally distributed with isotropic covariance with respect to the Riemannian metric of the space. Then, from observations, we aim to infer parameters of this metric. In the Large Deformation Diffeomorphic Metric Mapping (LDDMM) setting, this can be framed as inferring parameters of the kernel mappingbetween the dual Lie algebra and the Lie algebra of the diffeomorphism group that acts on the domain containing the landmarks. We achieve this by deriving a scheme for Monte Carlo simulation of Brownian landmark bridges conditioned on hitting the observed landmark configurations. Based on the Euclidean diffusion bridge simulation scheme of Delyon and Hu , we can compute expectation over bridges using the correction factor of a guided diffusion process to obtain the transition density of the Brownian motion. From this, we can take derivatives to obtain an iterative optimization algorithm for the most likely parameters. The scheme applies to the situation when the landmark configurations are considered observed at a fixed positive time . The time interval will generally be sufficiently large that many time discretization points are needed to accurately represent the stochastic process.
We begin in Section 2 with a short survey of LDDMM landmark geometry, metric estimation, large deformation stochastics, and uses of Brownian motion in shape analysis. In Section 3, we derive a scheme for simulating Brownian landmark bridges. We apply this scheme in Section 4 to derive an inference algorithm for estimating parameters of the metric. Numerical examples are presented in Section 5 before the paper ends with concluding remarks.
2 Landmarks Manifolds and Stochastic Landmark Dynamics
We start with a short survey of landmark geometry with the LDDMM framework as derived in papers including [24, 7, 11, 5]. The framework applies to general shape spaces though we focus on configurations of landmarks . We denote the resulting manifold . Two sets of shapes are in LDDMM matched by minimizing the energy functional
The parameter of
is a time-dependent vector fieldthat via a reconstruction equation
generates a corresponding time-dependent flow of diffeomorphisms . The endpoint diffeomorphism move the landmarks through the action of on . The right-most term of (1) measures the dissimilarity between and weighted by a factor . In the landmark case, the squared Euclidean distance when considering the landmarks elements of is often used here.
The Lagrangian on is often in the form with the -pairing and being a differential operator. Because can formally be considered the Lie algebra of , puts the dual Lie algebra into correspondence with by the mapping , . The inverse of the mapping arise from the Green’s function of , written as the kernel . Such defines a right-invariant inner product on the tangent bundle that descends to a Riemannian metric on . Because can be considered a subset of using the representation above, the metric structure can be written directly as a cometric
using the kernel evaluated on for two covectors
. The kernel is often specified directly in the form
kernels . One choice of is the Gaussian kernel with matrix specifying the spatial
correlation structure, and a scaling of the general kernel amplitude.
Estimating parameters of , with as above and the entries of or , has to our knowledge previously only been treated for landmarks in the small-deformation setting . While a linear vector space probability distribution is mapped to the manifold with small deformations, this paper concerns the situation when the probability distribution is constructed directly from the Riemannian metric on the nonlinear space . The approach has similarities with the estimation procedures derived in  where a metric on a finite dimensional Lie group is estimated to optimize likelihood of data on a homogeneous space arising as the quotient of the group by a closed subgroup. Though the landmark space can be represented as with the landmark isotropy subgroup , the approach of  can not directly be applied because of the infinite dimensionality of .
2.1 Brownian Motion
A diffusion processes on a Riemannian manifold is said to be a Brownian motion if its generator is with being the Laplace-Beltrami operator of the metric . Such processes can be constructed in several ways, see e.g. . By isometrically embedding in a Euclidean space , the process can be constructed as a process in that will stay on a.s. The process can equivalently be characterized in coordinates as being solution to the Itô integral
is a square root of the cometric tensor, and the drift term arise from contraction of the Christoffel symbols with the cometric. The noise term is infinitesimal increments of an -valued Brownian motion . Equivalently, the Brownian motion can be constructed as a hypoelliptic diffusion processes in the orthonormal frame bundle where a set of globally defined horizontal vector fields gives the Stratonovich integral equation
Note the sum over the horizontal fields . The process where is the bundle map is then a Brownian motion. This is known as the Eells-Elworthy-Malliavin construction of Brownian motion. The fields evaluated at model infinitesimal parallel transport of the vectors comprising the frame in the direction of the th frame vector, see e.g. .
While Brownian motion is per definition isotropic with equal variation in all directions, data with nontrivial covariance can be modeled by defining the SDE (5) in the larger frame bundle [17, 21] using nonorthonormal frames to model the square root of the local covariance structure. In this setup, the inference problem consists of finding the starting point of the diffusion and the square root covariance matrix. Estimators are defined via a Frechét mean like minimization in with square distances used as proxy for the negative log-transition density. In this paper, we remove this proxy by approximating the actual transition density, but only in the isotropic Brownian motion case.
2.2 Large Deformation Stochastics
Several papers have recently derived models for Brownian motion  and stochastic dynamics in shape analysis and for landmark manifolds. [25, 26] considered stochastic shape evolution by adding finite and infinite dimensional noise in the momentum equation of the dynamics. In , noise is added to the momentum equation to make the landmark dynamics correspond to a type of heat bath appearing in statistical physics. In [3, 4] a stochastic model for shape evolution is derived that descends to the landmark space in the same fashion as the right-invariant LDDMM metric descends to . The fundamental structure is here the momentum map that is preserved by the introduction of right-invariant noise. The approach is linked to parametric SDEs in fluid dynamics  and stochastic coadjoint motion .
3 Brownian Bridge Simulation
Brownian motion can be numerically simulated on using the coordinate Itô form (4). With a standard Euler discretization, the scheme becomes
with time discretization , and discrete noise , . Alternatively, a Heun scheme for discrete integration of the Stratonovich equation (5) results in
Because the horizontal fields represent infinitesimal parallel transport, they can be expressed using the Christoffel symbols of . The Christoffel symbols for the landmark metric are derived in  from which they can be directly implemented or implicitly retrieved from an automatic symbolic differentiation as done in the experiments in Section 5.
3.1 Bridge Sampling
The transition density of a Brownian motion evaluated at at time can be informally obtained by taking an expectation to get the “mass” of those of the sample paths hitting at time . We write for the process conditioned on hitting a.s. at . Computing the expectation analytically is in nonlinear situations generally intractable. Instead, we wish to employ a Monte Carlo approach and thus derive a method for simulating from . For this, we employ the bridge sampling scheme of Delyon and Hu . We first describe the framework for a general diffusion process in Euclidean space before using it directly on the landmark manifold .
be an valued Itô diffusion with invertible diffusion field . In order to sample from the conditioned process , , a modified processes is in  constructed by adding an extra drift term to the process giving the new process
As , the attraction term becomes increasingly strong forcing the processes to hit at a.s. It can be shown that the process exists when , and are with bounded derivatives. The process is then absolutely continuous with respect to the conditioned process . The Radon-Nikodym derivative between the laws and is
with denoting expectation with respect to , and the correction factor defined as the limit of
Here , , and quadratic variation is denoted by . Then and
for . The fact that the diffusion field must be invertible for the scheme to work as outlined here can be seen explicitly from the use of the inverse of and in these equations.
We can use the guided process (9) to take conditional expectation for general measurable functions on the Wiener space of paths by sampling from . Taking the particular choice of the constant function, the expression
for the transition density of arise as shown in . Note that both the leading factors and the correction factor are dependent on the diffusion field , the starting point of the diffusion, and the drift . Again, we can approximate the expectation in (11) by sampling from .
3.2 Landmark Bridge Simulation
Because the landmark manifold has a global chart on from the standard representation of each landmark position in , we can conveniently apply the bridge construction of . Writing the Itô coordinate form of the Brownian motion (4) in the form (8), we have and giving the guided SDE
The attraction term is the difference between the current landmark configuration and the target configuration . The transition density becomes
where we use the subscript to emphasize the dependence on the parameters and the kernel . As above, the expectation can be approximated by drawing samples from and evaluating .
A similar scheme is used for the bridge simulation of the stochastic coadjoint processes of [3, 4]. In these cases, the flow is hypoelliptic in the phase space and observations are partial in that only the landmark positions are observed. The momenta are unobserved. In addition, the fact that the landmarks can carry a large initial momentum necessitates a more general form of the guidance term that takes into account the expected value of of the process at time given the current time position and momentum of the process.
4 Inference Algorithm
Given a set of i.i.d. observations of landmark configurations, we assume the configurations are distributed according to the time transition distribution of a Brownian motion on started at . We now intend to infer parameters of the model. With the metric structure on given by (3) and kernel of the form , , parameters are the starting position , , and , i.e. .
The likelihood of the model given the data with respect to the Lebesgue measure on is
Using our ability to approximate (13) by bridge sampling, we aim to find a maximum-likelihood estimate (MLE) . We do this by a gradient based optimization on , see Algorithm 1. Note that the likelihood and thus the MLE of are dependent on the chosen background measure, in this case coming from the canonical chart on .
The inference Algorithm 1 optimizes the likelihoodand the observation time are considered missing data. The -step of the EM algorithm would then involve the expectation of the landmark bridges conditioned on the data with formally denoting a likelihood of an unconditioned sample path . This approach is used in e.g. . While natural to formulate, the approach involves the likelihood of a stochastic path which is only defined for finite time discretizations. In addition, the expected correction factor that arise when using the guided process in the estimation appears as a normalization factor in the EM -function. This can potentially make the scheme sensitive to the stochasticity in the Monte Carlo sampling of the expected correction . While the differences between these approaches needs further investigation, we hypothesize that direct optimization of the likelihood is superior in the present context.
Instead of taking expectations over , we can identify the most probable path of the conditioned process . This results in the Onsager-Machlup functional . In , a different definition is given that, in the isotropic Brownian motion situation, makes the set of Riemannian geodesics from to equal to the set of most probable paths of the conditioned process . The sample Frechét mean
is in that case formally also a minimizer of the negative log-probability of the most probable path to the data. Given that we are now able to approximate the density function of the Brownian motion, the MLE of the likelihood (14) with respect to is equivalent to
Compared to (15), the negative log-probability of the data is here minimized instead of the squared geodesic distance. The estimator (16) can therefore be considered a transition density equivalent of the sample Frechét mean.
5 Numerical Experiments
We here present examples of the method on simulated landmark configurations, and an application of the method and algorithm to landmarks annotated on cardiac images of left ventricles. We aim for testing the ability of the algorithm to infer the parameters of the model given samples. We here take the first steps in this direction and leave a more extensive simulation study to future work.
For the simulated data, we compare the results against the true values used in the simulation. In addition, we do simple model checking for both experiments by simulating with the estimated parameters and comparing the per-landmark sample mean and covariance.
We use code based on the Theano library for the implementation, in particular the symbolic expression and automatic derivative facilities of Theano. The code used for the experiments is available in the software package Theano Geometry http://bitbucket.org/stefansommer/theanogeometry. The implementation and the use of Theano for differential geometry applications including landmark manifolds is described in .
With 10 landmarks arranged in an ellipse configuration , we sample 64 samples from the transition distribution at time of a Brownian motion started at , see Figure 2. Parameters for the kernel are with set to the average inter-point distance in , and the amplitude parameter .
We run Algorithm 1 with initial conditions for the pointwise mean of the samples. The parameter evolution trough the iterative optimization and the result of the inference can be seen in Figure 3. The algorithm is able estimate the initial configuration and the parameters of and with a reasonable precision. The sample per-landmark covariance as measured on a new set of simulated data with the estimated parameters is comparable to the per-landmark covariance of the original dataset.
5.1 Left Cardiac Ventricles
To exemplify the approach on real data, we here use a set of landmarks obtained from annotations of the left ventricle in 14 cardiac images . Each ventricle is annotated with sets of landmarks from which we select 17 from each configuration for use in this experiment. Figure 4 shows an annotated image along with the sets of annotations for all images.
Figure 5 shows the results of the inference algorithm with setup equivalent to Figure 3. While the parameters converges during the iterative optimization, we here have no ground-truth comparison. A subsequent sampling using the estimated parameters allows comparison of the per-landmark sample covariance. While the new sample covariance in magnitude and to some degree shape corresponds to the sample covariance from the original data, the fact that the Brownian motion is isotropic forces the covariance to be equivalent for all landmarks as measured by the Riemannian landmark metric. Including anisotropic covariance in the distribution or the right-invariant stochastics of [3, 4] would allow the per-landmark covariance to vary and result in a closer fit.
In the paper, we have derived a method for maximum likelihood estimation of parameters for the starting point of landmark Brownian motions and for the Riemannian metric structure specified from the kernel . Using the guided process scheme of  for sampling conditioned Brownian bridges, the transition density is approximated by Monte Carlo sampling. With this approximation of the data likelihood, we use a gradient based iterative scheme to optimize parameters. We show on synthetic and real data sets the ability of the method to infer the underlying parameters of the data distribution and hence the metric structure of the landmark manifold.
A direct extension of the method presented here is to generalize to the anisotropic normal distributions  defined via Brownian motions in the frame bundle . This would allow differences in the per-landmark covariance and thus improve results on datasets such as the presented left ventricle annotations. Due to the hypoellipticity of the anisotropic flows that must be conditioned on hitting fibers in above points , further work is necessary to adapt the scheme presented here to the anisotropic case.
Acknowledgements We are grateful for the use of the cardiac ventricle dataset provided by Jens Chr. Nilsson and Bjørn A. Grønning, Danish Research Centre for Magnetic Resonance (DRCMR).
-  Allassonnière, S., Amit, Y., Trouvé, A.: Towards a coherent statistical framework for dense deformable template estimation. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69(1), 3–29 (Feb 2007)
-  Arnaudon, A., Castro, A.L., Holm, D.D.: Noise and dissipation on coadjoint orbits. arXiv:1601.02249 [math-ph, physics:nlin] (Jan 2016)
-  Arnaudon, A., Holm, D.D., Pai, A., Sommer, S.: A Stochastic Large Deformation Model for Computational Anatomy. Information Processing in Medical Imaging (IPMI) (2017)
-  Arnaudon, A., Holm, D.D., Sommer, S.: A Geometric Framework for Stochastic Shape Analysis. submitted, arXiv:1703.09971 [cs, math] (Mar 2017)
-  Beg, M.F., Miller, M.I., Trouvé, A., Younes, L.: Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. IJCV 61(2) (2005)
-  Delyon, B., Hu, Y.: Simulation of conditioned diffusion and application to parameter estimation. Stochastic Processes and their Appl. 116(11), 1660–1675 (2006)
-  Dupuis, P., Grenander, U., Miller, M.I.: Variational Problems on Flows of Diffeomorphisms for Image Matching (1998)
-  Fujita, T., Kotani, S.i.: The Onsager-Machlup function for diffusion processes. Journal of Mathematics of Kyoto University 22(1), 115–130 (1982)
-  Holm, D.D.: Variational principles for stochastic fluid dynamics. Proceedings. Mathematical, Physical, and Eng. Sciences / The Royal Society 471(2176) (2015)
-  Hsu, E.P.: Stochastic Analysis on Manifolds. American Mathematical Soc. (2002)
-  Joshi, S., Miller, M.: Landmark matching via large deformation diffeomorphisms. Image Processing, IEEE Transactions on 9(8), 1357–1370 (2000)
Kühnel, L., Arnaudon, A., Sommer, S.: Deep Learning Numerics for Differential Geometry. in preparation (2017)
-  Markussen, B.: Large deformation diffeomorphisms with application to optic flow. Comput. Vis. Image Underst. 106(1), 97–105 (Apr 2007)
-  Marsland, S., Shardlow, T.: Langevin equations for landmark image registration with uncertainty. arXiv:1605.09276 [math] (May 2016), arXiv: 1605.09276
-  Micheli, M.: The differential geometry of landmark shape manifolds: metrics, geodesics, and curvature. Ph.D. thesis, Brown University, Providence, USA (2008)
-  Papaspiliopoulos, O., Roberts, G.O.: Importance sampling techniques for estimation of diffusion models. In: Statistical Methods for Stochastic Differential Equations. Chapman & Hall/CRC Press (2012)
-  Sommer, S.: Anisotropic Distributions on Manifolds: Template Estimation and Most Probable Paths. In: Information Processing in Medical Imaging. Lecture Notes in Computer Science, vol. 9123, pp. 193–204. Springer (2015)
-  Sommer, S.: Anisotropically Weighted and Nonholonomically Constrained Evolutions on Manifolds. Entropy 18(12), 425 (Nov 2016)
-  Sommer, S., Jacobs, H.O.: Reduction by Lie Group Symmetries in Diffeomorphic Image Registration and Deformation Modelling. Symmetry 7(2), 599–624 (2015)
-  Sommer, S., Joshi, S.: Brownian Bridge Simulation and Metric Estimation on Lie Groups and Homogeneous Spaces. in preparation (2017)
-  Sommer, S., Svane, A.M.: Modelling Anisotropic Covariance using Stochastic Development and Sub-Riemannian Frame Bundle Geometry. Accepted for publication in Journal of Geometric Mechanics, arXiv:1512.08544 [math, stat] (2016)
-  Stegmann, M.B., Fisker, R., Ersbøll, B.K.: Extending and applying active appearance models for automated, high precision segmentation in different image modalities. Proc. 12th - SCIA 2001, Bergen, Norway pp. 90–97 (2001)
-  Team, T.T.D.: Theano: A Python framework for fast computation of mathematical expressions. arXiv:1605.02688 [cs] (May 2016), arXiv: 1605.02688
Trouvé, A.: An Infinite Dimensional Group Approach for Physics Based Models in Patterns Recognition (1995)
-  Trouvé, A., Vialard, F.X.: Shape splines and stochastic shape evolutions: A second order point of view. Quarterly of Applied Mathematics 70(2), 219–251 (2012)
-  Vialard, F.X.: Extension to infinite dimensions of a stochastic second-order model associated with shape splines. Stochastic Processes and their Applications 123(6), 2110–2157 (Jun 2013)