Colin J. Cotter

Diffeomorphic 3D Image Registration via Geodesic Shooting Using an Efficient Adjoint Calculation

In the context of large deformations by diffeomorphisms, we propose a new diffeomorphic registration algorithm for 3D images that performs the optimization directly on the set of geodesic flows. The key contribution of this work is to provide an accurate estimation of the so-called initial momentum, which is a scalar function encoding the optimal deformation between two images through the Hamiltonian equations of geodesics. Since the initial momentum has proven to be a key tool for statistics on shape spaces, our algorithm enables more reliable statistical comparisons for 3D images.Our proposed algorithm is a gradient descent on the initial momentum, where the gradient is calculated using standard methods from optimal control theory. To improve the numerical efficiency of the gradient computation, we have developed an integral formulation of the adjoint equations associated with the geodesic equations.We then apply it successfully to the registration of 2D phantom images and 3D cerebral images. By comparing our algorithm to the standard approach of Beg et al. (Int. J. Comput. Vis. 61:139–157, 2005), we show that it provides a more reliable estimation of the initial momentum for the optimal path. In addition to promising statistical applications, we finally discuss different perspectives opened by this work, in particular in the new field of longitudinal analysis of biomedical images.

A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling

18.02.14 KB. Ok to add accepted version to spiral, Elsevier says ok while mandate not enforced.

Mixed finite elements for numerical weather prediction

We show how mixed finite element methods that satisfy the conditions of finite element exterior calculus can be used for the horizontal discretisation of dynamical cores for numerical weather prediction on pseudo-uniform grids. This family of mixed finite element methods can be thought of in the numerical weather prediction context as a generalisation of the popular polygonal C-grid finite difference methods. There are a few major advantages: the mixed finite element methods do not require an orthogonal grid, and they allow a degree of flexibility that can be exploited to ensure an appropriate ratio between the velocity and pressure degrees of freedom so as to avoid spurious mode branches in the numerical dispersion relation. These methods preserve several properties of the C-grid method when applied to linear barotropic wave propagation, namely: (a) energy conservation, (b) mass conservation, (c) no spurious pressure modes, and (d) steady geostrophic modes on the f-plane. We explain how these properties are preserved, and describe two examples that can be used on pseudo-uniform grids: the recently-developed modified RTk-Q(k-1) element pairs on quadrilaterals and the BDFM1- P 1 DG element pair on triangles. All of these mixed finite element methods have an exact 2:1 ratio of velocity degrees of freedom to pressure degrees of freedom. Finally we illustrate the properties with some numerical examples.

Multisymplectic formulation of fluid dynamics using the inverse map

We construct multisymplectic formulations of fluid dynamics using the inverse of the Lagrangian path map. This inverse map, the ‘back-to-labels’ map, gives the initial Lagrangian label of the fluid particle that currently occupies each Eulerian position. Explicitly enforcing the condition that the fluid particles carry their labels with the flow in Hamiltons principle leads to our multisymplectic formulation. We use the multisymplectic one-form to obtain conservation laws for energy, momentum and an infinite set of conservation laws arising from the particle relabelling symmetry and leading to Kelvins circulation theorem. We discuss how multisymplectic numerical integrators naturally arise in this approach.

Continuous and Discrete Clebsch Variational Principles

Abstract The Clebsch method provides a unifying approach for deriving variational principles for continuous and discrete dynamical systems where elements of a vector space are used to control dynamics on the cotangent bundle of a Lie group via a velocity map. This paper proves a reduction theorem which states that the canonical variables on the Lie group can be eliminated, if and only if the velocity map is a Lie algebra action, thereby producing the Euler–Poincaré (EP) equation for the vector space variables. In this case, the map from the canonical variables on the Lie group to the vector space is the standard momentum map defined using the diamond operator. We apply the Clebsch method in examples of the rotating rigid body and the incompressible Euler equations. Along the way, we explain how singular solutions of the EP equation for the diffeomorphism group (EPDiff) arise as momentum maps in the Clebsch approach. In the case of finite-dimensional Lie groups, the Clebsch variational principle is discretized to produce a variational integrator for the dynamical system. We obtain a discrete map from which the variables on the cotangent bundle of a Lie group may be eliminated to produce a discrete EP equation for elements of the vector space. We give an integrator for the rotating rigid body as an example. We also briefly discuss how to discretize infinite-dimensional Clebsch systems, so as to produce conservative numerical methods for fluid dynamics.

Computational Modes and Grid Imprinting on Five Quasi-Uniform Spherical C Grids

AbstractCurrently, most operational forecasting models use latitude–longitude grids, whose convergence of meridians toward the poles limits parallel scaling. Quasi-uniform grids might avoid this limitation. Thuburn et al. and Ringler et al. have developed a method for arbitrarily structured, orthogonal C grids called TRiSK, which has many of the desirable properties of the C grid on latitude–longitude grids but which works on a variety of quasi-uniform grids. Here, five quasi-uniform, orthogonal grids of the sphere are investigated using TRiSK to solve the shallow-water equations.Some of the advantages and disadvantages of the hexagonal and triangular icosahedra, a “Voronoi-ized” cubed sphere, a Voronoi-ized skipped latitude–longitude grid, and a grid of kites in comparison to a full latitude–longitude grid are demonstrated. It is shown that the hexagonal icosahedron gives the most accurate results (for least computational cost). All of the grids suffer from spurious computational modes; this is especiall...

A FRAMEWORK FOR MIMETIC DISCRETIZATION OF THE ROTATING SHALLOW-WATER EQUATIONS ON ARBITRARY POLYGONAL GRIDS ∗

Accurate simulation of atmospheric flow in weather and climate prediction models requires the discretization of the governing equations to have a number of desirable properties. Although these properties can be achieved relatively straightforwardly on a latitude-longitude grid, they are much more challenging on the various quasi-uniform spherical grids that are now under consideration. A recently developed scheme—called TRiSK—has these desirable properties on grids that have an orthogonal dual. The present work extends the TRiSK scheme into a more general framework suitable for grids that have a nonorthogonal dual, such as the equiangular cubed sphere. We also show that this framework fits within the wider framework of mimetic discretizations and discrete exterior calculus. One key ingredient is the definition of certain mapping operators that are discrete analogues of the Hodge star operator, enabling the definition of a compatible inner product. Discrete Coriolis terms are also included within the mimet...

The variational particle-mesh method for matching curves

Diffeomorphic matching (only one of several names for this technique) is a technique for non-rigid registration of curves and surfaces in which the curve or surface is embedded in the flow of a time series of vector fields. One seeks the flow between two topologically-equivalent curves or surfaces which minimizes some metric defined on the vector fields, i.e. the flow closest to the identity in some sense. In this paper, we describe a new particle-mesh discretization for the evolution of the geodesic flow and the embedded shape. Particle-mesh algorithms are very natural for this problem because Lagrangian particles (particles moving with the flow) can represent the movement of the shape whereas the vector field is Eulerian and hence best represented on a static mesh. We explain the derivation of the method, and prove the conservation properties: the discrete method has a set of conserved momenta corresponding to the particle-relabelling symmetry which converge to conserved quantities in the continuous problem. We also introduce a new discretization for the geometric current matching condition of Vaillant and Glaunes (2005 Surface matching via currents Proc. Conf. IPMI pp 381–92). We illustrate the method and the derived properties with numerical examples.

The remapped particle-mesh semi-Lagrangian advection scheme

We describe the remapped particle-mesh advection method, a new mass-conserving method for solving the density equation which is suitable for combining with semi-Lagrangian methods for compressible flow applied to numerical weather prediction. In addition to the conservation property, the remapped particle-mesh method is computationally efficient and at least as accurate as current semi-Lagrangian methods based on cubic interpolation. We provide results of tests of the method in the plane, results from incorporating the advection method into a semi-Lagrangian method for the rotating shallow-water equations in planar geometry, and results from extending the method to the surface of a sphere.

A primal-dual mimetic finite element scheme for the rotating shallow water equations on polygonal spherical meshes

A new numerical method is presented for solving the shallow water equations on a rotating sphere using quasi-uniform polygonal meshes. The method uses special families of finite element function spaces to mimic key mathematical properties of the continuous equations and thereby capture several desirable physical properties related to balance and conservation. The method relies on two novel features. The first is the use of compound finite elements to provide suitable finite element spaces on general polygonal meshes. The second is the use of dual finite element spaces on the dual of the original mesh, along with suitably defined discrete Hodge star operators to map between the primal and dual meshes, enabling the use of a finite volume scheme on the dual mesh to compute potential vorticity fluxes. The resulting method has the same mimetic properties as a finite volume method presented previously, but is more accurate on a number of standard test cases.