Geometry and Shape Spaces

Contents

• Introduction and overview
• 0. Review: gradients, flows and reproducing kernels
• 1. Riemannian geometry and Kendall’s shape space
• 2. Lie groups, actions and metrics induced by deformations
• 3. Groups of diffeomorphisms built from kernels
• 4. Landmark matching: Hamiltonian dynamics and geodesic shooting
• 5. Matching without correspondences: measures, currents and varifolds
• 6. Images, binary regions and tensor fields
• 7. Metamorphoses: deforming and changing intensity
• 8. Statistics on shape spaces and anatomical variability
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with an M1 in mathematics, computer science or physics, at ease with differential calculus and linear algebra, who wants to compare, average and model shapes such as anatomical structures, curves, surfaces and images.

Empirical first: add a large smooth displacement to a grid and it folds, with a Jacobian of -1.06 in Code 3.4, while the flow of the same field never does; average two elongated diffusion tensors entry by entry and their volume grows by a factor 3.6 (Code 2.6). Linear operations fail on shapes, and this course builds the geometry that replaces them.

The course derives distances between shapes from the cost of the deformations that carry one shape to another. The induced Riemannian metrics give geodesics, exponential maps, registration algorithms (ch. 4–7) and statistical shape models (ch. 8). Related pages: Deformable Models and Geodesic Methods, Kernel Methods, Geometric Data Analysis and Medical Image Analysis.

What you will be able to do

• Define a Riemannian metric, compute geodesics, exponential and logarithm maps, and explain where normal coordinates stop being valid.
• Compute shape distances in Kendall’s space and between SPD matrices, and explain why entrywise averages of shapes or tensors fail.
• Build the metric induced on objects by a group of deformations, and derive its cometric and its Hamiltonian geodesic equations.
• Construct a group of diffeomorphisms from a reproducing kernel space of vector fields, and prove that its flows never fold space.
• Implement geodesic shooting and inexact matching for landmarks, and check energy conservation and the endpoint condition numerically.
• Compare point clouds, curves and surfaces without correspondences with kernel norms on measures, currents and varifolds.
• Derive and implement the LDDMM gradient for images, and explain why image momenta live on edges.
• Choose an action and a reorientation strategy for vector and tensor fields.
• Formulate image and landmark metamorphoses, and predict the effect of the intensity weight \sigma.
• Estimate a template and its modes of variation from a population, and state the limits of these estimates.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com sets the objective of an introductory exposition of a geometric viewpoint on deformable models. The viewpoint consists in building Riemannian metrics on spaces of objects, regarded as manifolds of finite or infinite dimension, from the study of the natural infinitesimal deformations acting on these objects.

The objects listed are n-tuples of points, curves, binary regions, surfaces, tensor fields, grey-level images and measures.

The themes (“Thèmes abordés”) say that building and studying these metrics, computing geodesics in shape spaces, and understanding local exponential maps matter both in theory and in practice. They name the concrete questions that this viewpoint unifies, announce medical illustrations on anatomical variability, and stress the new tools and algorithms for image processing. Each item is matched below.

The Université Paris-Saclay course detail gives the volume (21 h of lectures and 9 h of labs) and the nature of the assessment, but no programme items.

Beyond the syllabus, this page adds a review chapter (ch. 0), a research chapter for self-study (ch. 9), runnable code for every core chapter, and a written-exam proposal.

Where this course fits in the programme

Official prerequisites: master-mva.com states that the necessary mathematical tools, notably Riemannian geometry, Lie groups and Hamiltonian mechanics, are introduced along the way; no other prerequisite is listed. The instructors listed on master-mva.com are Alain Trouvé (Centre Borelli, ENS Paris-Saclay) and Joan Alexis Glaunès (Université Paris Cité). The course carries the label of the Santé track.

• Refreshers — Linear Algebra, Analysis and Optimization for Hilbert spaces, the spectral theorem and gradients; Probability and Statistics for the Gaussian models of ch. 8; Scientific Python, NumPy and PyTorch for the code.
• Semester 1 courses that prepare it — Geometric Data Analysis (kernel distances, shapes as measures), Optimal Transport for Machine Learning (the other family of data terms), Geometry Processing and Geometric Deep Learning (meshes), Medical Image Analysis (registration and computational anatomy) and Stochastic Calculus in Machine Learning (for the stochastic models of ch. 9).
• Related courses of semester 2 — Deformable Models and Geodesic Methods for Image Analysis, Kernel Methods for Machine Learning, Point Clouds and 3D Modelling, Deep Learning for Medical Imaging and Representation Learning for Computer Vision and Medical Imaging.
• Afterwards — the Research Internship; the Hospital Immersion of the Santé track gives a clinical context for the medical applications.

Course outline

• Review: gradients, flows and reproducing kernels — gradients and inner products, Euler–Lagrange and Hamilton’s equations, flows as diffeomorphisms, minimal-norm interpolation in an RKHS.
• Riemannian geometry and Kendall’s shape space — metrics, geodesics, the exponential map, quotients by isometries, and planar shape space.
• Lie groups, actions and metrics induced by deformations — invariant metrics (SPD matrices), projected metrics, cometrics, and the Euler–Poincaré equation.
• Groups of diffeomorphisms built from kernels — admissible spaces, flows, Trouvé’s group G_V, variations of flows, EPDiff.
• Landmark matching — reduction to a Hamiltonian system, geodesic shooting, the exponential chart, inexact matching.
• Matching without correspondences — point clouds as measures, curves and surfaces as currents and varifolds.
• Images, binary regions and tensor fields — actions, normal momenta, the LDDMM gradient, tensor reorientation.
• Metamorphoses — deformation plus intensity change, for images and landmarks.
• Statistics on shape spaces — Fréchet means, tangent PCA, generative deformable templates, medical illustrations.
• Research directions and open problems — metrics on curves, curvature, modular and stochastic deformations, new data terms, computation and learning (self-study).

Session plan

30 h in class = 10 sessions of 3 h.

The volume is the Université Paris-Saclay programme table: 21 h of lectures and 9 h of labs. master-mva.com gives no session organisation. The split into ten sessions is a proposal: session 1 is a lecture, and each later session ends with about one hour of lab on the code of the chapter.

Level 1 = foundations, 2 = core methods, 3 = advanced material close to research; the level never decreases from one session to the next.

Self-study, outside the 30 h (about 4 h in total): ch. 9 Research directions and open problems — metrics on curves, curvature and statistics, modular and stochastic deformations, data terms, computation, learning (≈ 4 h, with Ex. 16).

How to study this module

• Read ch. 0 even if it looks familiar: Proposition 0.2 (gradients depend on the inner product) and Proposition 0.10 (minimal-norm interpolation) are used throughout ch. 3–8.
• Follow one thread through the page: the cometric K(q) of Theorem 2.8 becomes the kernel matrix of landmarks (ch. 4), ch. 5 keeps the same shooting and only changes the data term, and metamorphosis only adds \sigma^2 I to the cometric (ch. 7).
• After each session, run the code of the chapter and change one parameter (the kernel width, the noise level, the weight \sigma) to see what breaks; Exercise 14 does this systematically.
• Do the exercises listed in the session plan before the next session; the ★★★ exercises are good starting points for a project.
• Python environment: NumPy and PyTorch on CPU. Every code block runs in under 30 seconds on a laptop with synthetic data; GPU kernel libraries are only needed for the larger projects.

0. Review: gradients, flows and reproducing kernels

The whole course builds distances between shapes out of the cost of deforming them. Three M1 tools carry that construction: gradients, which depend on the chosen inner product; ordinary differential equations, whose flows are diffeomorphisms; and Hilbert spaces of functions with a reproducing kernel. This chapter states them in the exact form used later, with the proofs that later chapters rely on.

Gradients depend on the inner product

The differential is intrinsic, but the gradient is not: it changes when the inner product changes. This simple fact explains why every algorithm of the course smooths its descent directions.

In ch. 3–7, K is a smoothing kernel. The V-gradient is therefore the L^2-gradient smoothed by K, and descent along it keeps the deformations smooth. Numerically, it acts as a preconditioner.

Calculus of variations and Hamiltonian mechanics

This is all the Hamiltonian mechanics the course needs. In ch. 4 the cometric K(q) will be a kernel matrix, and Proposition 0.5 becomes the landmark matching algorithm.

Flows of vector fields

The Jacobian stays positive however large the velocity is. This is the root of the large deformation models of the syllabus: the flow of a smooth field never folds space (Code 3.4).

Hilbert spaces and reproducing kernels

The main example is the Gaussian kernel k(x, y) = \exp(-|x - y|^2 / \sigma^2). Its Fourier transform is a positive Gaussian, so it is positive definite (Bochner’s theorem), and even strictly: its Gram matrices at distinct points are invertible.

1. Riemannian geometry and Kendall’s shape space

A space of shapes needs three things: a notion of length for a continuous change of shape, the shortest paths it defines, and a way to linearize the space around a reference shape. Riemannian geometry provides all three. This chapter sets them up in finite dimension and builds the first shape space, Kendall’s space of landmark configurations modulo translation, scale and rotation.

Manifolds, metrics and geodesics

All finite-dimensional spaces of this course are of this kind. Spaces of diffeomorphisms, curves and images are infinite-dimensional; their geometry is built by hand in ch. 3–7.

Minimizing the energy rather than the length is the standard choice: the energy is smooth and convex in the velocity, and the parametrization comes out right automatically.

Geodesics are the critical points of the energy. They minimize length between points that are close enough, but not necessarily between distant ones: on the sphere, a great circle stops minimizing after half a turn.

The exponential map

The coordinates v = \log_x(y) are the normal coordinates at x. They are the linear coordinates in which shape statistics are computed (ch. 8). Two limits matter in practice.

• The exponential map may stop being injective: beyond the cut locus, geodesics from x no longer minimize.
• Its differential may become singular: such points are conjugate to x. On the unit sphere, the antipode of x is both, at distance \pi.

This theorem fails in infinite dimension (Atkin, 1975). The existence of optimal deformations in ch. 3, 4 and 7 is therefore proved directly, by compactness arguments.

Quotients by a group of isometries

A shape is what remains of a configuration once position, size and orientation are forgotten. Mathematically, this is a quotient by a group that acts by isometries.

When a Lie group acts freely and properly by isometries on a Riemannian manifold, the quotient is a manifold. Its metric is the restriction of the metric to horizontal vectors, those orthogonal to the orbit. Horizontal geodesics project onto geodesics of the quotient (a Riemannian submersion; see Lee, 2018).

Kendall (1984) showed that \Sigma_2^k is the complex projective space \mathbb{C}P^{k-2} with the Fubini–Study metric. Its diameter is \pi / 2.

For triangles (k = 3) it is a sphere of radius 1/2. The two mirror-image equilateral triangles are its poles, and the collinear triangles form its equator.

After the optimal rotation, \langle z, w \rangle is real and positive, and the great circle from z to the rotated w is horizontal. Its projection is therefore a geodesic of shape space.

import numpy as np

def preshape(X):
    """Centre and scale k planar landmarks; return them as a unit complex vector."""
    z = X[:, 0] + 1j * X[:, 1]
    z = z - z.mean()
    return z / np.linalg.norm(z)

def kendall_dist(X, Y):
    """Distance in Kendall's shape space of planar configurations (Theorem 1.11)."""
    z, w = preshape(X), preshape(Y)
    return np.arccos(min(1.0, abs(np.vdot(z, w))))

def geodesic(X, Y, t):
    """Point at time t of the shape geodesic: rotate Y onto X, then follow the great circle."""
    z, w = preshape(X), preshape(Y)
    c = np.vdot(w, z)
    w = w * c / abs(c)                       # optimal rotation of w (Procrustes)
    rho = np.arccos(min(1.0, abs(np.vdot(z, w))))
    g = (np.sin((1 - t) * rho) * z + np.sin(t * rho) * w) / np.sin(rho)
    return np.c_[g.real, g.imag]

equi = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3) / 2]])
flat = np.array([[0, 0], [1, 0], [2, 0]])
right = np.array([[0, 0], [1, 0], [0, 1]])
mirror = equi * [1, -1]                      # the reflected equilateral triangle

print(f"d(equilateral, collinear)  = {kendall_dist(equi, flat):.4f}   (pi/4 = {np.pi/4:.4f})")
print(f"d(equilateral, right)      = {kendall_dist(equi, right):.4f}")
print(f"d(equilateral, mirror)     = {kendall_dist(equi, mirror):.4f}   (pi/2 = {np.pi/2:.4f})")

th = 0.7                                     # similarity invariance
Rot = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
moved = 3.0 * right @ Rot.T + [5.0, -2.0]
print(f"d(equilateral, moved right) = {kendall_dist(equi, moved):.4f}")

mid = geodesic(equi, flat, 0.5)              # the midpoint is halfway on both sides
print(f"midpoint: d to each end = {kendall_dist(equi, mid):.4f}, {kendall_dist(mid, flat):.4f}")

In three dimensions the rotation group no longer acts freely: a collinear configuration is fixed by the rotations about its line, and the quotient has singular points. The optimal rotation is then given by the singular value decomposition (Procrustes analysis; Gower, 1975).

2. Lie groups, actions and metrics induced by deformations

Kendall’s space came from a group of symmetries: moves that do not change the shape. The course’s main construction uses a group differently. The group deforms the objects, the cost of an infinitesimal deformation is fixed once and for all, and the metric on objects is the cheapest way to produce a given change. This chapter carries out the construction in finite dimension, where each step can be checked by hand.

Matrix Lie groups

The velocity \xi = \dot g g^{-1} is the right-trivialized, or Eulerian, velocity. It is the finite-dimensional model of the flow equation \partial_t \varphi_t = v_t \circ \varphi_t of ch. 3, where v_t is the velocity field seen at the current position of the points.

Group actions

When the group is a symmetry: invariant metrics

A metric on Q is invariant if every map q \mapsto g \cdot q is an isometry. Kendall’s metric is invariant under rotations. The second classical example, used for diffusion tensor fields in ch. 6, lives on symmetric positive definite (SPD) matrices.

import numpy as np

def spd_fun(S, f):
    """Apply f to the eigenvalues of a symmetric matrix S."""
    w, U = np.linalg.eigh(S)
    return (U * f(w)) @ U.T

def dist_ai(A, B):
    """Affine-invariant distance ||log(A^{-1/2} B A^{-1/2})||_F (Proposition 2.5)."""
    Ai = spd_fun(A, lambda w: w ** -0.5)
    return np.linalg.norm(spd_fun(Ai @ B @ Ai, np.log))

def geo_ai(A, B, t):
    """Affine-invariant geodesic A^{1/2} (A^{-1/2} B A^{-1/2})^t A^{1/2}."""
    Ah, Ai = spd_fun(A, np.sqrt), spd_fun(A, lambda w: w ** -0.5)
    return Ah @ spd_fun(Ai @ B @ Ai, lambda w: w ** t) @ Ah

def rot(deg):
    a = np.radians(deg)
    return np.array([[np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]])

A = np.diag([4.0, 0.25])                     # an elongated "diffusion tensor"
B = rot(60) @ A @ rot(60).T                  # the same tensor turned by 60 degrees

mid_euc = (A + B) / 2
mid_ai = geo_ai(A, B, 0.5)
mid_le = spd_fun((spd_fun(A, np.log) + spd_fun(B, np.log)) / 2, np.exp)
for name, M in [("Euclidean", mid_euc), ("affine-invariant", mid_ai), ("log-Euclidean", mid_le)]:
    print(f"{name:17s} midpoint: det = {np.linalg.det(M):.3f}, eigenvalues = {np.round(np.linalg.eigvalsh(M), 3)}")

rng = np.random.default_rng(0)
G = rng.standard_normal((2, 2))              # any invertible matrix acts by S -> G S G^T
print(f"d(A, B) = {dist_ai(A, B):.4f}, d(GAG^T, GBG^T) = {dist_ai(G @ A @ G.T, G @ B @ G.T):.4f}")
print(f"d(A, mid) = {dist_ai(A, mid_ai):.4f}, d(mid, B) = {dist_ai(mid_ai, B):.4f}")

The log-Euclidean mean, \exp\big(\tfrac12 (\log A + \log B)\big), is a cheaper alternative that also avoids swelling (Arsigny, Fillard, Pennec and Ayache, 2007).

When the group deforms: projected metrics

Choose an inner product \langle \cdot, \cdot \rangle on the Lie algebra and declare the cost of a path g(t) to be \int |\dot g g^{-1}|^2 dt. This defines a right-invariant metric on G: right multiplication by a fixed element does not change the Eulerian velocity. The metric on objects follows.

To move n points independently, the group must be infinite-dimensional, and the linear kernel must be replaced by a kernel whose Gram matrices are invertible, such as the Gaussian. This is the group of diffeomorphisms of ch. 3, and its cometric on landmarks is the Gaussian kernel matrix (ch. 4).

The Euler–Poincaré equation

Arnold (1966) showed that the Euler equations of a perfect incompressible fluid are of this form, for the group of volume-preserving diffeomorphisms with the L^2 metric. On the group of diffeomorphisms with the metric of V, the equation is called EPDiff (ch. 3). Its conservation form says that the whole geodesic is determined by its initial momentum. This is the basis of geodesic shooting.

3. Groups of diffeomorphisms built from kernels

To move any configuration of points, any curve or any image, the deforming group must be infinite-dimensional: it is a group of diffeomorphisms of \mathbb{R}^d. Its “Lie algebra” is a Hilbert space V of smooth vector fields, and the norm of V is the cost of an infinitesimal deformation. The construction (Trouvé, 1998; Dupuis, Grenander and Miller, 1998) follows ch. 2 step by step, but each finite-dimensional fact has to be proved again.

Admissible spaces of vector fields

The standard choice is K(x, y) = k(x, y)\, I_d with the Gaussian k(x, y) = \exp(-|x - y|^2 / \sigma^2). Sobolev spaces H^s with s > d/2 + 1 are also admissible, by the Sobolev embedding theorem.

In operator form, \|v\|_V^2 = \langle L v, v \rangle_{L^2}, with L the inverse of the convolution by K. For Sobolev spaces, L = (I - \alpha \Delta)^s.

The momentum of v is m = L v, so that v = K m. In ch. 4 the momentum is a sum of Dirac masses, and in ch. 6 it is concentrated on the edges of an image.

Flows and the group G_V

The simplest alternative to a flow is the “small deformation” model x \mapsto x + v(x). It is guaranteed to be a diffeomorphism when \|Dv\|_\infty < 1, and it folds as soon as \det(I + Dv) changes sign. The next code compares the two for a single Gaussian push.

import numpy as np

sigma, s = 0.25, 0.6                         # kernel width and strength of the push

def v(x):
    """Stationary field s * k(x, 0) * (1, 0): one Gaussian kernel at the origin, pushing right."""
    k = np.exp(-(x ** 2).sum(-1, keepdims=True) / sigma ** 2)
    return s * k * np.array([1.0, 0.0])

n = 201
g = np.linspace(-1, 1, n)
X = np.stack(np.meshgrid(g, g, indexing="ij"), -1)   # grid of points, shape (n, n, 2)
h = g[1] - g[0]

def min_jacobian(Phi):
    """Smallest determinant of the finite-difference Jacobian of a map sampled on the grid."""
    dx, dy = np.gradient(Phi, h, axis=(0, 1))
    return (dx[..., 0] * dy[..., 1] - dx[..., 1] * dy[..., 0]).min()

linear = X + v(X)                            # "small deformation" model x -> x + v(x)

Phi, steps = X.copy(), 100                   # flow of v for unit time (Euler steps)
for _ in range(steps):
    Phi = Phi + v(Phi) / steps

print(f"bound 1 - s*sqrt(2)/sigma*exp(-1/2) = {1 - s * np.sqrt(2) / sigma * np.exp(-0.5):.3f}")
print(f"min det D(x + v(x)) = {min_jacobian(linear):.3f}")
print(f"min det D(phi_1)    = {min_jacobian(Phi):.3f}")
print(f"displacement of the origin: linear {linear[100, 100, 0]:.3f}, flow {Phi[100, 100, 0]:.3f}")

G_V is not a Lie group in the classical sense, because composition is not smooth for the natural topology. Yet the calculus of ch. 2 goes through, with one tool: the derivative of the flow with respect to the velocity.

Variations of a flow

Every gradient computed in ch. 4–7 is an application of this formula, followed by Proposition 0.2 to turn an L^2-gradient into a V-gradient.

Geodesics: EPDiff and momentum conservation

For a diffeomorphism \varphi and a vector field w, write \operatorname{Ad}_\varphi w = (D\varphi\, w) \circ \varphi^{-1} , the field w transported by \varphi.

This is Theorem 2.10 for the group G_V, where \operatorname{Ad}_\varphi replaces \eta \mapsto g \eta g^{-1}. A complete proof requires care with the infinite-dimensional setting (Miller, Trouvé and Younes, 2006; Younes, 2010). For landmarks, ch. 4 checks it by a direct computation (Proposition 4.4).

Holm, Marsden and Ratiu (1998) derived the same equation as the Euler–Poincaré equation of the diffeomorphism group.

4. Landmark matching: Hamiltonian dynamics and geodesic shooting

The simplest objects that G_V deforms are n-tuples of points, the landmarks that an anatomist places on recognizable features. Matching two sets of labelled landmarks (Joshi and Miller, 2000) reduces the infinite-dimensional problem of ch. 3 to a Hamiltonian system in \mathbb{R}^{2nd}, which can be solved exactly and fast. The curve, surface and point-cloud methods of ch. 5 reuse this reduction unchanged.

From diffeomorphisms to a finite-dimensional problem

Landmark matching is therefore a geodesic problem on the finite-dimensional manifold Q_n, with metric G(q) = K(q)^{-1}. This is exactly the projected metric of Theorem 2.8, with A_q v = v \cdot q and A_q^* p = \sum_j K(\cdot, q_j) p_j.

Geodesic shooting and the exponential map

Given an initial momentum p_0, integrating Theorem 4.3 on [0, 1] is geodesic shooting. It computes the exponential map q(1) = \exp_{q_0}\big( K(q_0)\, p_0 \big).

By Proposition 1.6, p_0 \mapsto q(1) is a local diffeomorphism near 0, with differential K(q_0): the momentum is a linear chart of landmark space around the template. Every other point z of the ambient space is carried along by \dot z = \sum_j K(z, q_j)\, p_j. This is the dense matching from landmarks of the syllabus.

For large momenta the chart can fail: the differential of the shooting map may become singular (conjugate points). Landmark space has sectional curvature of both signs (Micheli, Michor and Mumford, 2012), so no global chart should be expected.

Inexact matching

Measured landmarks are noisy, and exact matching would follow the noise. Inexact matching trades the deformation cost against the residual. Since optimal paths are geodesics, it suffices to optimize over the initial momentum: J(p_0) = H(q_0, p_0) + \frac{1}{2\lambda^2}\, |q(1) - y|^2, \qquad q(1) = \exp_{q_0}\big( K(q_0)\, p_0 \big).

The gradient of J with respect to p_0 is obtained by integrating the linearized Hamiltonian system backwards in time from this endpoint condition (the adjoint method). Reverse-mode automatic differentiation through the integrator computes the same quantity, and is what the next code uses.

import math
import torch

torch.set_default_dtype(torch.float64)
torch.set_num_threads(1)
sigma, lam, steps = 0.25, 0.1, 20            # kernel width, data noise level, time steps

def kmat(x, y):
    return torch.exp(-((x[:, None, :] - y[None, :, :]) ** 2).sum(-1) / sigma ** 2)

def hamiltonian(q, p):
    return 0.5 * ((p @ p.T) * kmat(q, q)).sum()

def field(q, p):
    """Right-hand side of Hamilton's equations for the Gaussian kernel (Theorem 4.3)."""
    K = kmat(q, q)
    dq = K @ p
    w = K * (p @ p.T) * (2 / sigma ** 2)                 # w_ij = (2/sigma^2) k_ij <p_i, p_j>
    dp = w.sum(1, keepdim=True) * q - w @ q              # sum_j w_ij (q_i - q_j)
    return dq, dp

def shoot(q, p, extra=None):
    """Geodesic shooting with the midpoint rule; extra points are carried by the flow."""
    dt = 1.0 / steps
    for _ in range(steps):
        dq, dp = field(q, p)
        qm, pm = q + 0.5 * dt * dq, p + 0.5 * dt * dp
        if extra is not None:
            em = extra + 0.5 * dt * kmat(extra, q) @ p
            extra = extra + dt * kmat(em, qm) @ pm
        dq, dp = field(qm, pm)
        q, p = q + dt * dq, p + dt * dp
    return q, p, extra

n = 20
th = torch.linspace(0, 2 * math.pi, n + 1)[:-1]
q0 = 0.5 * torch.stack([th.cos(), th.sin()], 1)                       # template: a circle
r = 0.5 + 0.15 * torch.cos(3 * th)
y = torch.stack([r * th.cos(), r * th.sin()], 1) + torch.tensor([0.1, 0.0])   # labelled target

p0 = torch.zeros_like(q0, requires_grad=True)
opt = torch.optim.LBFGS([p0], max_iter=100, tolerance_grad=1e-10, line_search_fn="strong_wolfe")
def closure():
    opt.zero_grad()
    q1, _, _ = shoot(q0, p0)
    J = hamiltonian(q0, p0) + ((q1 - y) ** 2).sum() / (2 * lam ** 2)   # inexact matching
    J.backward()
    return J
opt.step(closure)

with torch.no_grad():
    g = torch.linspace(-1, 1, 41)
    grid = torch.stack(torch.meshgrid(g, g, indexing="ij"), -1).reshape(-1, 2)
    q1, p1, phi = shoot(q0, p0, grid)
    H0, H1 = hamiltonian(q0, p0).item(), hamiltonian(q1, p1).item()
    phi, h = phi.reshape(41, 41, 2), (g[1] - g[0]).item()
    dx = (phi[2:, 1:-1] - phi[:-2, 1:-1]) / (2 * h)
    dy = (phi[1:-1, 2:] - phi[1:-1, :-2]) / (2 * h)
    det = dx[..., 0] * dy[..., 1] - dx[..., 1] * dy[..., 0]
    res = q1 - y
    print(f"RMS landmark gap: {((q0 - y) ** 2).sum(1).mean().sqrt():.4f} -> {(res ** 2).sum(1).mean().sqrt():.4f}")
    print(f"geodesic length sqrt(2H) = {math.sqrt(2 * H0):.4f}; relative drift of H = {abs(H1 - H0) / H0:.1e}")
    K1 = kmat(q1, q1)                                  # endpoint condition, read on velocities
    print(f"endpoint condition: |K(p(1) + (q(1) - y)/lam^2)| / |K p(1)| = {(K1 @ (p1 + res / lam ** 2)).norm() / (K1 @ p1).norm():.1e}")
    print(f"Jacobian determinant of phi_1 on the grid: min {det.min():.3f}, max {det.max():.3f}")

The endpoint condition is checked on velocities, K(q_1)\, p(1), rather than on momenta. The Gram matrix of nearby landmarks is badly conditioned, and momentum components along its small eigenvectors move nothing; the optimizer leaves them unresolved.

Two parameters must be chosen with care. The noise level \lambda is a distance, in the units of the landmarks. The width \sigma is the scale of coherent motion (Example 4.5); anatomical data often call for a sum of Gaussians of several widths (Risser, Vialard, Wolz, Murgasova, Holm and Rueckert, 2011).

5. Matching without correspondences: measures, currents and varifolds

Landmarks require point-to-point correspondences. Point clouds, curves and surfaces extracted from scans have different numbers of points and no labels. The remedy is to embed every object in a Hilbert space where a deformed template and a target can be compared directly. The Hamiltonian machinery of ch. 4 is then reused, with only the data term changed.

Point clouds as measures

This data term was introduced for unlabelled point sets by Glaunès, Trouvé and Younes (2004). In machine learning the same quantity is the maximum mean discrepancy between two samples.

Matching a template cloud \mu to a target \nu then minimizes J(p_0) = H(q_0, p_0) + \frac{1}{2\lambda^2}\, \| \varphi_* \mu - \nu \|_{W^*}^2 , where q_0 are the points of \mu. Shooting is that of ch. 4; only the gradient of the endpoint term changes. This is the dense matching from point clouds of the syllabus.

The weights a_i matter. With a_i = 1/n the measure is a probability, and two samplings of the same shape at different densities are close. With a_i = 1, a denser cloud has more mass and the distance sees it.

Curves and surfaces as currents

A measure forgets the geometry of a curve: it sees the points, not the tangents. A curve is better described by how vector fields flow through it.

Currents were introduced in this setting for surfaces by Vaillant and Glaunès (2005), and for curves by Glaunès, Qiu, Miller and Younes (2008). They solve the correspondence problem and are insensitive to the sampling. Their linearity, however, has a price.

Varifolds: forgetting the orientation

All terms of a varifold norm between nearby, nearly parallel elements are positive: there is no cancellation. In Example 5.5 the pair has the same varifold norm, 2(1 + k_W), whatever the orientations. Varifolds also handle meshes whose orientation is inconsistent or undefined.

Normal cycles add curvature information (Roussillon and Glaunès, 2016), and optimal transport gives another family of data terms (ch. 9).

import math
import torch

torch.set_default_dtype(torch.float64)
torch.set_num_threads(1)
sigma, tau, lam, steps = 0.25, 0.2, 0.05, 20   # deformation kernel, varifold kernel, noise, steps

def kmat(x, y, s):
    return torch.exp(-((x[:, None, :] - y[None, :, :]) ** 2).sum(-1) / s ** 2)

def field(q, p):
    K = kmat(q, q, sigma)
    w = K * (p @ p.T) * (2 / sigma ** 2)
    return K @ p, w.sum(1, keepdim=True) * q - w @ q

def shoot(q, p):
    dt = 1.0 / steps
    for _ in range(steps):
        dq, dp = field(q, p)
        dq2, dp2 = field(q + 0.5 * dt * dq, p + 0.5 * dt * dp)
        q, p = q + dt * dq2, p + dt * dp2
    return q, p

def segments(x):
    """Centres, unit tangents and lengths of the segments of a closed polygon."""
    t = torch.roll(x, -1, 0) - x
    length = t.norm(dim=1, keepdim=True)
    return (x + torch.roll(x, -1, 0)) / 2, t / length, length

def vf_inner(A, B):
    """Varifold inner product (Definition 5.6): sum l_i l_j k(c_i, c_j) <u_i, u_j>^2."""
    (c1, u1, l1), (c2, u2, l2) = A, B
    return ((l1 @ l2.T) * kmat(c1, c2, tau) * (u1 @ u2.T) ** 2).sum()

def vf_dist2(x, y):
    A, B = segments(x), segments(y)
    return vf_inner(A, A) - 2 * vf_inner(A, B) + vf_inner(B, B)

def closed_curve(n, r, shift=0.0):
    th = torch.linspace(0, 2 * math.pi, n + 1)[:-1] + shift
    return torch.stack([r(th) * th.cos(), r(th) * th.sin()], 1)

x0 = closed_curve(40, lambda t: 0.5 + 0 * t)                            # template: circle, 40 points
y = closed_curve(64, lambda t: 0.45 + 0.2 * torch.cos(2 * t), shift=1.0)  # target: 64 points, no labels
y = y * torch.tensor([1.0, 0.6]) + torch.tensor([0.1, 0.05])             # flattened into a peanut-like shape

p0 = torch.zeros_like(x0, requires_grad=True)
opt = torch.optim.LBFGS([p0], max_iter=60, line_search_fn="strong_wolfe")
def closure():
    opt.zero_grad()
    x1, _ = shoot(x0, p0)
    J = 0.5 * ((p0 @ p0.T) * kmat(x0, x0, sigma)).sum() + vf_dist2(x1, y) / (2 * lam ** 2)
    J.backward()
    return J
opt.step(closure)

with torch.no_grad():
    x1, _ = shoot(x0, p0)
    def mean_gap(a, b):                                # symmetric mean distance to the nearest vertex
        D = ((a[:, None] - b[None]) ** 2).sum(-1).sqrt()
        return 0.5 * (D.min(1).values.mean() + D.min(0).values.mean())
    print(f"varifold distance^2: {vf_dist2(x0, y):.4f} -> {vf_dist2(x1, y):.5f}")
    print(f"mean gap to the nearest vertex of the other curve: {mean_gap(x0, y):.4f} -> {mean_gap(x1, y):.4f}")
    print(f"deformation cost sqrt(2H) = {math.sqrt(((p0 @ p0.T) * kmat(x0, x0, sigma)).sum()):.4f}")

The cost of these data terms is a double sum over the elements of the two shapes. With 10^5 triangles it is computed on a GPU without storing the kernel matrix, with libraries such as KeOps (ch. 9).

6. Images, binary regions and tensor fields

For grey-level images the object is a function, the group acts by composition, and many velocity fields change nothing at all: those that slide along the level lines. The structure of ch. 2 still applies. It yields the large-deformation image registration algorithm of Beg, Miller, Trouvé and Younes (2005), and explains why its momenta live on the edges of the image.

Actions of diffeomorphisms on fields

Each formula is an action. For images, \varphi \cdot (\psi \cdot I) equals I \circ (\varphi \circ \psi)^{-1}; for vector and tensor fields, the chain rule D(\varphi \circ \psi) = (D\varphi \circ \psi)\, D\psi does the same job.

For diffusion tensor images, the full action also stretches the tensors. Alexander, Pierpaoli, Basser and Gee (2001) proposed instead to apply only the rotation part R = (D\varphi D\varphi^\top)^{-1/2} D\varphi of the Jacobian (finite strain), or to follow the deformed principal direction (preservation of principal direction).

Momenta of images are thus concentrated where the gradient is large: on the edges. For a binary region, \nabla \mathbf{1}_A = -n\, \sigma_{\partial A}, the outward normal times the surface measure of the boundary. The momentum is then normal to the boundary.

To match a region A to a region B, the natural data term is the area of the symmetric difference: \| \mathbf{1}_{\varphi(A)} - \mathbf{1}_B \|_{L^2}^2 = | \varphi(A) \,\triangle\, B | .

Write A_t = \varphi_t(A), whose boundary moves with velocity v. The area of the symmetric difference can be rewritten as | A_t \,\triangle\, B | = |B| + \int_{A_t} \big( 1 - 2\, \mathbf{1}_B \big)\, dx .

When \partial A_t meets \partial B transversally, the transport formula for a moving domain then gives \frac{d}{dt}\, | A_t \,\triangle\, B | = \int_{\partial A_t} \big( 1 - 2\, \mathbf{1}_B \big)\, (v \cdot n)\, d\sigma .

The L^2-gradient is therefore a normal vector measure carried by the boundary. By Proposition 0.2, the V-gradient is this measure smoothed by K: \nabla_V\, | A_t \,\triangle\, B | = K\big( (1 - 2\, \mathbf{1}_B)\, n\, \sigma_{\partial A_t} \big) .

A descent step pushes the boundary outwards where it lies inside B, and inwards where it lies outside B. Matching regions thus amounts to evolving their boundary curves or surfaces with normal momenta.

Check on a disk A of radius r inside a larger concentric disk B, under the dilation v(x) = x. At t = 0 both sides equal -2\pi r^2.

Large deformation matching of images

At a critical point, v_t = K(P_t\, \nabla J^0_t), with the scalar momentum P_t = \frac{1}{\lambda^2}\, |\det D\varphi_{t,1}|\, \big( J^0_t - J^1_t \big).

As Proposition 6.3 predicts, the momentum is a scalar times the image gradient. It is transported as a density, \partial_t P + \operatorname{div}(P v) = 0, while the image is transported, \partial_t J^0 + \nabla J^0 \cdot v = 0.

The pair (I_0, P_0) therefore determines the whole geodesic: this is geodesic shooting for images (Miller, Trouvé and Younes, 2006). Each subject of a population is then summarized by a scalar image P_0 on the template, the representation used in ch. 8.

import math
import torch
import torch.nn.functional as F

torch.set_default_dtype(torch.float64)
torch.set_num_threads(1)
N, T, sig_v, sig_d = 48, 8, 0.15, 0.05       # grid size, time steps, kernel width, data noise
h = 2 / (N - 1)                               # pixel size on [-1, 1]^2
g = torch.linspace(-1, 1, N)
Y, X = torch.meshgrid(g, g, indexing="ij")
ident = torch.stack([X, Y], -1)[None]         # grid_sample convention: (x, y) in [-1, 1]

def soft(sd):                                 # smoothed indicator of {sd < 0}
    return torch.sigmoid(-sd / 0.04)

I0 = soft(torch.sqrt(X ** 2 + Y ** 2) - 0.45)                    # a disk
d1 = torch.sqrt((X - 0.3) ** 2 + Y ** 2) - 0.32
d2 = torch.sqrt((X + 0.3) ** 2 + Y ** 2) - 0.32
I1 = soft(torch.minimum(d1, d2))                                 # a peanut (union of two disks)

f = torch.fft.fftfreq(N, d=h) * 2 * math.pi
Kh = torch.exp(-(f[:, None] ** 2 + f[None, :] ** 2) * sig_v ** 2 / 4)   # Gaussian kernel, Fourier side

def smooth(u):                                # v = K * u, channel by channel
    return torch.fft.ifft2(torch.fft.fft2(u) * Kh).real

def warp(img, grid):                          # img o grid, bilinear
    return F.grid_sample(img[None, None], grid, align_corners=True, padding_mode="border")[0, 0]

def flow_inverse(u):
    """phi_1^{-1} by the semi-Lagrangian rule psi_{t+dt} = psi_t o (id - dt v_t)."""
    psi = ident.clone()
    for t in range(T):
        v = torch.stack([smooth(u[t, 0]), smooth(u[t, 1])], -1)[None]
        disp = (psi - ident).permute(0, 3, 1, 2)
        back = ident - v / T
        psi = back + F.grid_sample(disp, back, align_corners=True, padding_mode="border").permute(0, 2, 3, 1)
    return psi

def energy(u):
    reg = sum((u[t] * torch.stack([smooth(u[t, 0]), smooth(u[t, 1])])).sum() for t in range(T)) * h ** 2 / T
    J = warp(I0, flow_inverse(u))
    return 0.5 * reg + ((J - I1) ** 2).sum() * h ** 2 / (2 * sig_d ** 2)

u = torch.zeros(T, 2, N, N, requires_grad=True)   # momenta: v_t = K u_t
opt = torch.optim.LBFGS([u], max_iter=40, line_search_fn="strong_wolfe")
def closure():
    opt.zero_grad()
    E = energy(u)
    E.backward()
    return E
opt.step(closure)

with torch.no_grad():
    psi = flow_inverse(u)
    J = warp(I0, psi)
    dx = (psi[0, 1:-1, 2:] - psi[0, 1:-1, :-2]) / (2 * h)       # derivative along x (columns)
    dy = (psi[0, 2:, 1:-1] - psi[0, :-2, 1:-1]) / (2 * h)       # derivative along y (rows)
    det = dx[..., 0] * dy[..., 1] - dx[..., 1] * dy[..., 0]
    print(f"relative L2 mismatch: {((I0 - I1).norm() / (I0 - I1).norm()):.3f} -> {((J - I1).norm() / (I0 - I1).norm()):.3f}")
    print(f"Jacobian determinant of phi_1^-1: min {det.min():.3f}, max {det.max():.3f}")

The largest Jacobian of \varphi_1^{-1}, about 10, sits at the waist of the peanut: the centre of the disk is squeezed there by that factor. The smallest, 0.33, lies inside the lobes, which the disk fills by expansion. The Gaussian convolution by FFT treats the image as periodic, which is harmless here because the shapes stay away from the border.

Tensor fields

For a diffusion tensor image \Sigma_0 and a target \Sigma_1, the data term compares the deformed field with the target, point by point. It can use the Frobenius norm or the affine-invariant distance of Proposition 2.5: \frac{1}{2\lambda^2} \int_\Omega d\big( (\varphi_1 \cdot \Sigma_0)(x),\, \Sigma_1(x) \big)^2\, dx .

The gradient is computed as in Theorem 6.5, with the chain rule through D\varphi_1 in the action. LDDMM for vector fields was developed by Cao, Miller, Winslow and Younes (2005), with an application to cardiac fibre orientation.

7. Metamorphoses: deforming and changing intensity

Diffeomorphisms only move what already exists. A template without a lesion can never be deformed into an image that has one, and the matching of ch. 6 then distorts the surrounding anatomy to compensate. Metamorphosis (Miller and Younes, 2001; Trouvé and Younes, 2005) lets the template itself change along the path, at a price. It defines a metric on the whole space of images, not only on the orbit of a template.

Definition and first properties

The residual Z is the material derivative: the change of intensity seen by a particle that moves with the deformation. A pure deformation, I_t = I_0 \circ \varphi_t^{-1}, has Z = 0. A pure change of intensity has v = 0 and Z = \partial_t I.

In Riemannian terms, the metric on images is |\delta I|_I^2 = \min_v \Big( \|v\|_V^2 + \frac{1}{\sigma^2}\, \|\delta I + \nabla I \cdot v\|_{L^2}^2 \Big) . It is the projected metric (Theorem 2.8) of the product of the group with the flat space of images.

Metamorphosis distances are therefore finite between any two images, unlike the LDDMM distance, which is infinite outside the orbit of the template. When \sigma \to \infty, intensity changes cost nothing and the optimal path tends to the cross-fade. When \sigma \to 0, they are forbidden and LDDMM is recovered.

The first equation is the optimality in v at fixed I. At time t the energy is a quadratic function of v, whose V-gradient is 2v + 2\sigma^{-2} K(Z \nabla I). The momentum is again a scalar times the image gradient (Proposition 6.3).

The second equation comes from the variations of I and says that the residual is transported as a density. Holm, Trouvé and Younes (2009) derived both equations from a general Euler–Poincaré theory of metamorphosis.

Landmark metamorphosis

For landmarks, metamorphosis lets each point drift with respect to the flow. The points obey \dot q_i = v_t(q_i) + \zeta_i, and the energy is \frac12 \int_0^1 \Big( \|v_t\|_V^2 + \frac{1}{\sigma^2}\, |\zeta|^2 \Big)\, dt .

The new cometric is invertible even when two landmarks coincide, so landmarks may now merge, and the space is completed. The shooting code of ch. 4 applies unchanged once \sigma^2 is added to the diagonal of the kernel matrix.

A discrete metamorphosis of signals

A path of images I_0, I_1, \dots, I_T with velocities v_0, \dots, v_{T-1} is discretized with a time step \delta t = 1/T. The residual uses the semi-Lagrangian transport of ch. 6: Z_t = \frac{I_{t+1} - I_t \circ (\mathrm{id} - \delta t\, v_t)}{\delta t} .

The intermediate images and the velocities are optimized together. This is the automatic metamorphosis between images of the syllabus.

import torch

torch.set_default_dtype(torch.float64)
torch.set_num_threads(1)
N, T, s_v = 128, 10, 0.08                     # samples, time steps, kernel width
x = torch.linspace(0, 1, N)
h, dt = 1 / (N - 1), 1 / T
K = torch.exp(-(x[:, None] - x[None, :]) ** 2 / s_v ** 2) * h   # v = K u (discrete kernel operator)

def bump(c, a=1.0, w=0.05):
    return a * torch.exp(-(x - c) ** 2 / (2 * w ** 2))

I0 = bump(0.4, w=0.06)                        # source: one bump
I1 = bump(0.52, w=0.06) + bump(0.85, 0.5, w=0.04)   # target: it has moved, a new one appeared

def resample(I, s):
    """Linear interpolation of the samples I at the positions s (clamped to [0, 1])."""
    s = s.clamp(0, 1) / h
    i = s.floor().clamp(max=N - 2).long()
    w = s - i
    return I[i] * (1 - w) + I[i + 1] * w

def costs(u, mid, sigma):
    """Discrete metamorphosis (Definition 7.1): deformation and intensity parts of the energy."""
    imgs = [I0] + list(mid) + [I1]
    e_def = e_int = 0.0
    for t in range(T):
        v = K @ u[t]
        Z = (imgs[t + 1] - resample(imgs[t], x - dt * v)) / dt   # residual intensity change
        e_def = e_def + 0.5 * dt * h * (u[t] @ v)
        e_int = e_int + 0.5 * dt * (Z ** 2).sum() * h / sigma ** 2
    return e_def, e_int, imgs

for sigma in (0.05, 0.5, 5.0):
    u = torch.zeros(T, N, requires_grad=True)
    mid = torch.stack([(1 - k * dt) * I0 + k * dt * I1 for k in range(1, T)]).requires_grad_(True)
    opt = torch.optim.LBFGS([u, mid], max_iter=300, line_search_fn="strong_wolfe")
    def closure():
        opt.zero_grad()
        e_def, e_int, _ = costs(u, mid, sigma)
        (e_def + e_int).backward()
        return e_def + e_int
    opt.step(closure)
    with torch.no_grad():
        e_def, e_int, imgs = costs(u, mid, sigma)
        track = " ".join(f"{x[im[:90].argmax()]:.2f}/{im[:90].max():.2f}" for im in imgs[::2])
        print(f"sigma = {sigma:4}: deformation {e_def:.3f}, intensity {e_int:.3f} | "
              f"main bump (position/height) at t = 0, 0.2, ..., 1: {track}")

The problem is not convex. When the two positions of a feature do not overlap at all, the cross-fade initialization gives no gradient towards transport, and the optimizer stays at the cross-fade. Coarse-to-fine schemes, which start with a wide kernel and smoothed images, are the standard remedy.

8. Statistics on shape spaces and anatomical variability

Computational anatomy asks how an organ varies across a population, and how a disease changes it. The constructions of ch. 1–7 provide the ingredients: a mean shape, or template; linear coordinates around it, given by the initial momenta of the exponential map; and probability models on these coordinates. This chapter assembles them into the probabilistic shape models announced by the syllabus.

Fréchet means

The Karcher mean is unique when the data lie in a geodesic ball that is small enough with respect to the curvature (Karcher, 1977; Afsari, 2011). In Kendall’s space this holds for populations whose shapes are not too different. In landmark spaces with curvature of both signs, uniqueness can fail for spread-out data (ch. 9).

Coordinates around the template

A principal mode w_k with variance \lambda_k is displayed by shooting from the template with \pm 2 \sqrt{\lambda_k}\, w_k. Principal geodesic analysis (Fletcher, Lu, Pizer and Joshi, 2004) replaces the linear subspaces of the tangent space by geodesic submanifolds; tangent PCA is its first-order approximation.

import math
import numpy as np
import torch

torch.set_default_dtype(torch.float64)
torch.set_num_threads(1)
sigma, steps, k, n_subj = 0.4, 20, 24, 40
rng = np.random.default_rng(0)

def kmat(x, y):
    return torch.exp(-((x[:, None, :] - y[None, :, :]) ** 2).sum(-1) / sigma ** 2)

def shoot(q, p):                                       # Hamiltonian geodesic shooting (ch. 4)
    dt = 1.0 / steps
    def field(q, p):
        K = kmat(q, q)
        w = K * (p @ p.T) * (2 / sigma ** 2)
        return K @ p, w.sum(1, keepdim=True) * q - w @ q
    for _ in range(steps):
        dq, dp = field(q, p)
        dq, dp = field(q + 0.5 * dt * dq, p + 0.5 * dt * dp)
        q, p = q + dt * dq, p + dt * dp
    return q

th = torch.linspace(0, 2 * math.pi, k + 1)[:-1]
q0 = torch.stack([0.6 * th.cos(), 0.25 * th.sin()], 1)          # template: an elongated outline
K0 = kmat(q0, q0)
mode1 = torch.stack([q0[:, 0], -q0[:, 1]], 1)                   # stretch along x, squeeze along y
mode2 = torch.stack([0 * q0[:, 0], q0[:, 0] ** 2 - 0.12], 1)    # bend: the ends go up
subjects, weights = [], []
for _ in range(n_subj):                                         # random initial velocities -> momenta
    a, b = rng.normal(0, 0.4), rng.normal(0, 0.25)
    weights.append((a, b))                                      # hidden from the analysis below
    p = torch.linalg.solve(K0, a * mode1 + b * mode2)          # p0 = K(q0)^{-1} u0
    subjects.append(shoot(q0, p).numpy() + 0.005 * rng.normal(size=(k, 2)))   # + observation noise

def pre(X):                                           # Kendall preshape, complex form (ch. 1)
    z = X[:, 0] + 1j * X[:, 1]
    z = z - z.mean()
    return z / np.linalg.norm(z)

def align(w, z):                                      # rotate w optimally onto z
    c = np.vdot(w, z)
    return w * c / abs(c)

def log_map(z, w):                                    # Riemannian log on the preshape sphere
    w = align(w, z)
    c = np.real(np.vdot(z, w))
    v = w - c * z
    nv = np.linalg.norm(v)
    return np.arccos(min(1.0, c)) * v / nv if nv > 1e-12 else 0 * z

def exp_map(z, v):
    t = np.linalg.norm(v)
    return np.cos(t) * z + np.sin(t) * v / t if t > 0 else z

Z = [pre(X) for X in subjects]
mu = Z[0]
for it in range(10):                                  # Fréchet mean by gradient descent (Prop. 8.2)
    step = np.mean([log_map(mu, w) for w in Z], axis=0)
    mu = exp_map(mu, step)
print(f"Frechet mean: norm of the mean log after 10 steps = {np.linalg.norm(step):.1e}")
print(f"distance from the Frechet mean to the template = {np.arccos(min(1, abs(np.vdot(mu, pre(q0.numpy()))))):.4f}")
V = np.array([log_map(mu, w) for w in Z])
V = np.c_[V.real, V.imag]                             # tangent vectors as real vectors
U, s, _ = np.linalg.svd(V - V.mean(0), full_matrices=False)
print("tangent PCA, share of variance of the first 4 modes:", np.round(s[:4] ** 2 / (s ** 2).sum(), 3))
scores = U[:, :2] * s[:2]                             # coordinates of the subjects on modes 1 and 2
corr = np.corrcoef(np.c_[scores, np.array(weights)].T)[:2, 2:]
print("|correlation| of modes 1, 2 with the hidden stretch, bend weights:", np.round(np.abs(np.diag(corr)), 3))

Probabilistic deformable templates

Under this prior, the initial velocities of the landmarks are the values at \bar q of a Gaussian random field with covariance kernel K. The whole field v_0 is the minimal-norm interpolant of these values (Proposition 0.10), that is, the conditional mean of that random field given its values at the landmarks. Minimizing the criterion alternately in the momenta (one inexact matching per subject, ch. 4) and in the template is the joint registration and template estimation of Joshi, Davis, Jomier and Gerig (2004).

Maximizing over the hidden momenta, instead of integrating them out, has a cost: the estimated template need not converge to the true one when the noise is not negligible. Allassonnière, Amit and Trouvé (2007) treat the deformations as random effects and estimate the template and the covariance by maximum likelihood, with consistency results. Allassonnière, Kuhn and Trouvé (2010) compute the estimate with a stochastic approximation EM algorithm.

Illustrations in medical imaging

These tools are used to analyse anatomical variability. Some representative studies follow.

• Csernansky and coauthors (1998) compared the hippocampi of patients with schizophrenia and of control subjects by high-dimensional brain mapping, and located the differences on the surface of the structure rather than in a single volume number.
• Durrleman, Pennec, Trouvé and Ayache (2009) built statistical models of sets of curves and surfaces with currents (ch. 5), with templates and modes of variation estimated from sulcal lines of the cortex.
• Longitudinal data, such as repeated scans of the same subject, are modelled by geodesic regression (Niethammer, Huang and Vialard, 2011; Fletcher, 2013) and by spatiotemporal atlases (Durrleman and coauthors, 2013).
• The open-source software Deformetrica (Bône, Louis, Martin and Durrleman, 2018) implements template estimation with landmarks, currents and varifolds.

The Medical Image Analysis course of semester 1 places these models among the other generative, geometric and biophysical approaches to medical images.

9. Research directions and open problems

This chapter is for self-study. It surveys active questions that extend the course, with entry points into the literature. Each section ends with problems that are still open or only partly solved.

Metrics directly on spaces of curves and surfaces

The course built metrics on curves from outside, through the deformations of the ambient space. One can also put a metric directly on parametrized curves, then quotient by reparametrizations. The simplest choice, the L^2 metric on velocities of curves, fails dramatically.

Michor and Mumford (2005, 2006) showed that the geodesic distance it induces on the space of unparametrized curves is identically zero: any two curves can be joined by paths of arbitrarily small length, which zigzag at high frequency. Sobolev metrics, which also penalize derivatives of the deformation along the curve, restore a positive distance.

A special first-order Sobolev metric, the elastic metric, becomes flat or spherical after a change of variables. Younes (1998) gave such a transform for plane curves. Srivastava, Klassen, Joshi and Jermyn (2011) introduced the square-root velocity transform, q = \dot c / \sqrt{|\dot c|}, under which the elastic metric becomes the L^2 metric.

Bauer, Bruveris and Michor (2014) survey these geometries and those of diffeomorphism groups. Open questions include the completeness of higher-order metrics on surfaces, the existence of optimal reparametrizations, and fast numerics for surfaces.

Curvature and the statistics of shapes

Micheli, Michor and Mumford (2012) computed the sectional curvature of landmark spaces with a kernel metric, through a formula in terms of the cometric. The curvature takes both signs. Positive curvature makes geodesics converge and Fréchet means non-unique for spread-out data.

The quotient structure creates a further bias. Devilliers, Allassonnière, Trouvé and Pennec (2017) showed that estimating a template by a Fréchet mean in a quotient space is not consistent when the noise is not negligible: the estimate moves away from the true template as the noise grows. How to correct this bias for realistic anatomical data remains open.

Constrained and modular deformations

The kernel of V encodes only smoothness. Biological growth and motion are better described by a few structured mechanisms: local scalings, rotations, growth along a direction. Gris, Durrleman and Trouvé (2018) introduced deformation modules, which build the velocity field from a small number of parametrized generators. The resulting geometry is sub-Riemannian.

Arguillère, Trélat, Trouvé and Younes (2015) recast shape deformation as optimal control, which covers constrained models and gives existence results and Pontryagin-type equations. Choosing the modules from data, rather than by hand, is an open problem.

Stochastic shape evolutions

Geodesics model smooth, deterministic evolutions; real growth is noisy. Trouvé and Vialard (2012) introduced shape splines and stochastic perturbations of the Hamiltonian dynamics of landmarks. Arnaudon, Holm and Sommer (2019) proposed a geometric framework where the noise enters the momentum equation in a way compatible with the group structure. The stochastic calculus course of semester 1 gives the tools.

Inference in these models is hard: the likelihood of observed shapes requires bridges of the stochastic dynamics, which are expensive to simulate. Estimating the noise structure from longitudinal data is an open problem.

Data attachment terms

Varifolds and currents are not the last word. Kaltenmark, Charlier and Charon (2017) unified currents and varifolds through oriented varifolds, whose kernel on directions interpolates between the two. Roussillon and Glaunès (2016) introduced kernel metrics on normal cycles, which see curvature and boundaries.

Feydy, Charlier, Vialard and Peyré (2017) used optimal transport as a data term for diffeomorphic registration. Sinkhorn divergences, introduced by Genevay, Peyré and Cuturi (2018), interpolate between optimal transport and kernel norms; Feydy and coauthors (2019) proved that they define positive, convex and smooth data terms; see the Optimal Transport for Machine Learning course. Which data term suits which anatomy, with which theoretical guarantees, is largely an empirical question today.

Computation at scale

Kernel sums dominate the cost: every shooting step and every data term is a sum over pairs of points. KeOps (Charlier, Feydy, Glaunès, Collin and Durif, 2021) computes such reductions on GPUs without storing the kernel matrix, with automatic differentiation. It makes LDDMM on meshes of 10^5 vertices practical.

Deformetrica (Bône, Louis, Martin and Durrleman, 2018) packages atlas estimation, and Geomstats (Miolane and coauthors, 2020) offers Riemannian tools, including Kendall shape spaces and SPD matrices, in Python.

Learning and shape spaces

Deep networks now predict registrations in one forward pass. Quicksilver (Yang, Kwitt, Styner and Niethammer, 2017) predicts the initial momentum of LDDMM from an image pair. VoxelMorph (Balakrishnan, Zhao, Sabuncu, Guttag and Dalca, 2019) learns a displacement field without supervision. Dalca, Balakrishnan, Guttag and Sabuncu (2019) obtain diffeomorphisms by integrating a predicted stationary velocity field.

These amortized methods are fast, but several questions remain open.

• Does the predicted map stay a diffeomorphism, and is the induced distance still a metric?
• How much does the network generalize across scanners, populations and pathologies?
• Can the kernel, or the deformation modules, be learned from data while keeping the geometric guarantees?

The courses on Deep Learning for Medical Imaging and Representation Learning for Computer Vision and Medical Imaging take up the learning side of these questions.

Open problems in brief

• Existence and regularity of metamorphosis geodesics for non-smooth images and for surfaces with topological changes.
• Consistent template estimation in quotient spaces, with explicit bias corrections.
• Fast, certified geodesic shooting for images in 3D at clinical resolution.
• Statistical tests on shape spaces with curvature, beyond the tangent-space linearization.
• Combining learned registration with the metric structure of G_V, so that learned maps define distances.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. The session plan says which ones to do after each session.

1. ★ Exponential and logarithm on the sphere. On the unit sphere, let x = (0, 0, 1). Compute \exp_x(v) for v = (\pi/3, 0, 0), then \log_x(y) for y = (0, 1, 0), and check that d(x, \exp_x(v)) = |v|.

   Solution: by Example 1.8, \exp_x(v) = \cos(\pi/3)\, x + \sin(\pi/3)\, e_1 = (0.866, 0, 0.5) . Its angle with x is \arccos 0.5 = \pi/3 = |v|.

   For y = (0, 1, 0), \theta = \pi/2 and y - \cos\theta\, x = y, so \log_x(y) = (0, \pi/2, 0).

2. ★ Kendall distance by hand. Compute the Kendall distance between the equilateral triangle (0, 1, e^{i\pi/3}) and the right isosceles triangle (0, 1, i).

   Solution: after centring, the equilateral triangle already has norm 1, and the right triangle w_c = \big( -(1+i), 2-i, -1+2i \big)/3 has norm 2/\sqrt3. The Hermitian product of the centred vectors is \langle z_c, w_c \rangle = \tfrac12 + \tfrac{\sqrt3}{3} + i\, \tfrac{\sqrt3}{6}, \qquad |\langle z_c, w_c \rangle| = 1.1154 . Dividing by 2/\sqrt3 gives 0.9659 = \cos(\pi/12), so the distance is \pi/12 = 0.2618, as in Code 1.12.

3. ★ Swelling. Let \Sigma_1 = I and \Sigma_2 = \operatorname{diag}(e^2, e^{-2}). Compute their affine-invariant distance, their Euclidean midpoint and its determinant, and their affine-invariant midpoint.

   Solution: the eigenvalues of \Sigma_1^{-1} \Sigma_2 are e^{\pm 2}, so d = \sqrt{4 + 4} = 2\sqrt2 = 2.828.

   The Euclidean midpoint is \operatorname{diag}(4.195, 0.568), with determinant 2.381, although both matrices have determinant 1. The affine-invariant midpoint is \operatorname{diag}(e, e^{-1}), with determinant 1.

4. ★ Moving together or apart. With the Gaussian kernel and two points at distance \sigma (\kappa = e^{-1}), compute the minimal RKHS norm of a function taking the values (1, 1), then (1, -1). Compare with Example 0.11.

   Solution: by Proposition 0.10, \|f^\star\|^2 = y^\top K^{-1} y. For y = (1, 1) it is 2/(1 + \kappa) = 1.462; for y = (1, -1) it is 2/(1 - \kappa) = 3.164.

   Coherent values are cheap and opposite values expensive; the case (1, 0) of Example 0.11, 1.156, lies in between.

5. ★ Infinitesimal actions. Compute the following infinitesimal actions (Definition 2.3 and the table that follows it).

   • 1. v \cdot I for the image I(x) = x_1^2 and the translation v = e_1.
   • 2. The action of the dilation v(x) = x on a density \rho of \mathbb{R}^d.
   • 3. \xi \cdot \Sigma for \Sigma = \operatorname{diag}(4, 1/4) and the rotation generator \xi with entries \xi_{12} = -1 and \xi_{21} = 1.

   Solution: (a) -\nabla I \cdot e_1 = -2 x_1. (b) -\mathrm{div}(\rho\, x) = -d\, \rho - x \cdot \nabla\rho.

   3. \xi\Sigma + \Sigma\xi^\top has zero diagonal and off-diagonal entries 4 - 1/4 = 3.75. A rotation changes an anisotropic tensor; for \Sigma = I the result would be \xi + \xi^\top = 0.

6. ★ Inexact matching of one landmark. One landmark at q_0, target y, kernel with k(x, x) = 1. Solve the inexact matching problem of ch. 4 in closed form and check Proposition 4.6.

   Solution: H = |p|^2/2, so p is constant and q(1) = q_0 + p_0. Minimizing \tfrac12 |p_0|^2 + |q_0 + p_0 - y|^2 / (2\lambda^2) gives p_0 = \frac{y - q_0}{1 + \lambda^2} . Then q(1) - y = -\lambda^2 (y - q_0)/(1 + \lambda^2), so -(q(1) - y)/\lambda^2 = p_0 = p(1).

7. ★ Two landmarks translated together. Two landmarks at distance r start with the same velocity u. Compute their initial momenta and the initial cost \|v_0\|_V^2. Is the momentum constant in time?

   Solution: K(q)\, p = (u, u) gives p_1 = p_2 = u/(1 + \kappa), with \kappa = e^{-r^2/\sigma^2}, and \|v_0\|_V^2 = 2|u|^2/(1 + \kappa).

   By Theorem 4.3, \dot p_1 = 2 \sigma^{-2} \kappa |p_1|^2 (q_1 - q_2) \neq 0: the momenta acquire components along q_1 - q_2, and the points drift slightly apart. The rigid translation is not a geodesic.

8. ★★ How much can a deformation compress? Using Proposition 0.7 and Definition 3.1, show that for v \in L^2([0, 1], V) \det D\varphi_1^v(x) \geq \exp\Big( -d\, c_V\, \Big( \int_0^1 \|v_t\|_V^2\, dt \Big)^{1/2} \Big) .

   Solution: the divergence is the trace of D v_t, so |\operatorname{div} v_t| = |\operatorname{tr} D v_t| \leq d\, \|D v_t\|_\infty \leq d\, c_V \|v_t\|_V . Insert this in the formula of Proposition 0.7, then bound \int_0^1 \|v_t\|_V\, dt by the Cauchy–Schwarz inequality.

   A deformation of length \ell cannot shrink volumes by more than a factor e^{-d c_V \ell}.

9. ★★ Noether’s theorem for landmarks. For the Gaussian kernel, show that the total momentum \sum_i p_i and, in the plane, the angular momentum \sum_i q_i \wedge p_i are conserved, where a \wedge b = a_1 b_2 - a_2 b_1.

   Solution: with w_{ij} = 2 \sigma^{-2} k(q_i, q_j)\, (p_i \cdot p_j), symmetric in (i, j), Theorem 4.3 reads \dot p_i = \sum_j w_{ij} (q_i - q_j). Summing over i, antisymmetric terms cancel, so \sum_i \dot p_i = 0.

   For the angular momentum, the derivative is \sum_{i,j} k(q_i, q_j)\, p_j \wedge p_i + \sum_{i,j} w_{ij}\, q_i \wedge (q_i - q_j) . The first sum is antisymmetric under the swap of i and j, and the second equals -\sum_{i,j} w_{ij}\, q_i \wedge q_j, also antisymmetric. Both vanish.

10. ★★ Head-on landmarks never collide. Let q_1 = -q_2 = (a, 0) and p_1 = -p_2 = (-b, 0) with a, b > 0. Show that the configuration stays symmetric, that H = b^2 (1 - \kappa) with \kappa = e^{-4a^2/\sigma^2}, and that a(t) decreases but never reaches 0.

    Solution: the symmetry (q, p) \mapsto (-q, -p) preserves the equations, so it holds for all times, with q_i and p_i on the first axis.

    Then \dot q_1 = (1 - \kappa)\, p_1, which gives \dot a = -(1 - \kappa)\, b and H = b^2 (1 - \kappa).

    Eliminating b with the conservation of H gives \dot a = -\sqrt{H (1 - \kappa)}. For small a, 1 - \kappa \approx 4a^2/\sigma^2, hence \dot a \approx -(2\sqrt{H}/\sigma)\, a: the gap decays exponentially and never closes in finite time.

11. ★★ EPDiff in one dimension. In dimension 1, Theorem 3.8 gives m_t = (\psi')^2\, m_0 \circ \psi with \psi = \varphi_t^{-1}. Derive the equation \partial_t m + v\, \partial_x m + 2\, (\partial_x v)\, m = 0.

    Solution: differentiating \varphi_t \circ \psi = \mathrm{id} in time gives \partial_t \psi = -v\, \psi', hence \partial_t \psi' = -v' \psi' - v\, \psi''. Then \partial_t m = 2 \psi' (\partial_t \psi')\, m_0 \circ \psi + (\psi')^2 (m_0' \circ \psi)\, \partial_t \psi = -2 v' m - v \big( 2 \psi' \psi''\, m_0 \circ \psi + (\psi')^3\, m_0' \circ \psi \big) . The bracket is exactly \partial_x m, so \partial_t m = -2 v' m - v\, \partial_x m. This is EPDiff with D v = \operatorname{div} v = v'.

12. ★★ The cross-fade bound for landmarks. For landmark metamorphosis (Proposition 7.4), show that the distance between two configurations x and y is at most |y - x| / \sigma. For a single landmark, compute it exactly.

    Solution: take v = 0 and \zeta = y - x constant: the energy is |y - x|^2 / (2\sigma^2), so d \leq |y - x| / \sigma.

    For one landmark, the cometric is the scalar 1 + \sigma^2; the geodesic is the straight line, and d = |y - x| / \sqrt{1 + \sigma^2}, which is indeed smaller than both |y - x| and |y - x| / \sigma.

13. ★★ The current of a small loop. Compute the squared current norm of a square of side a, discretized by its four sides, with a Gaussian kernel of width \tau. Show that it behaves like 4 a^4 / \tau^2 when a \ll \tau, and interpret.

    Solution: adjacent sides are orthogonal and contribute nothing. Each side contributes a^2 with itself, and each of the two pairs of opposite sides contributes twice -a^2 e^{-a^2/\tau^2}. The total is 4 a^2 \big( 1 - e^{-a^2/\tau^2} \big) \approx \frac{4 a^4}{\tau^2} . The norm, 2a^2/\tau, is proportional to the enclosed area: a closed loop much smaller than \tau is seen as a tiny vortex, not as a curve.

14. ★★ Choosing the kernel width (computational). Run Code 4.7 with \sigma = 0.1, 0.25 and 0.5. Report the geodesic length, the final gap and the minimal Jacobian, and explain the trends.

    Hints and expected result: the geodesic lengths are about 0.61, 0.48 and 0.78, and the minimal Jacobians about 0.14, 0.39 and 0.39. The final gaps stay near 0.001, except 0.005 for the widest kernel.

    With a narrow kernel, landmarks move almost independently: each pays for its own displacement, and the deformation concentrates near the points (smallest Jacobian). With a wide kernel, creating three lobes requires tearing neighbours apart, which Example 4.5 shows to be expensive. The intermediate width fits the scale of the shape.

15. ★★★ Atlas estimation (computational). Starting from the population of Code 8.4, estimate a template by alternating (a) an inexact LDDMM matching of the template to each subject (Code 4.7), and (b) an update of the template that moves it along the average momentum. Compare the estimated template with the true one and with the Kendall Fréchet mean, then perform PCA of the momenta with the inner product of Definition 8.3.

    Hints and expected result: at a stationary point the mean of the initial momenta, read as velocities, is close to zero, which is the analogue of Proposition 8.2. Report the Kendall distance of both estimates to the generating template, and how it changes with the observation noise; relate the trend to the bias discussed after Proposition 8.6. The first PCA mode should be the stretching mode.

16. ★★★ Currents versus varifolds on a fold (computational). Add to the target of Code 5.7 a thin spike: a narrow excursion that goes out and comes back along almost the same path. Match the circle with a current data term and with a varifold data term, and compare the results as the width of the spike decreases.

    Hints and expected result: the two sides of the spike have opposite tangents. Below the kernel width \tau they cancel as a current (Example 5.5), so current matching ignores the spike, while varifold matching grows it. Report the length of the matched curve inside the spike region against the spike width, for both data terms.

17. ★★★ Image metamorphosis in 2D. Extend Code 7.6 to two-dimensional images: a smooth disk as template, and the same disk, shifted, with a small bright spot added as target. Compare with LDDMM (Code 6.6) on the same pair.

    Hints and expected result: LDDMM cannot create the spot; it either ignores it or distorts the disk towards it, with large Jacobians. Metamorphosis transports the disk and creates the spot by intensity change. Plot the deformation and intensity parts of the energy against \sigma, and locate the value where the spot stops being produced by deformation.

Projects and assessment

The Université Paris-Saclay course detail describes the assessment as “Evaluation Continue non Intégrale” (partial continuous assessment). master-mva.com displays no assessment section, and no weights are published. The subjects below are proposals for the project part of such an assessment: each reproduces or extends a named research paper, or builds a complete pipeline on synthetic data with a measurable success criterion. The written exam is a self-assessment proposal, not an official modality.

Solo projects

1. [Solo · ≈ 20 h] Paper study: LDDMM for images. Reimplement the gradient descent of Beg, Miller, Trouvé and Younes (2005) with the explicit gradient of Theorem 6.5, not automatic differentiation.
   • Work: forward and backward flows on a 2D grid, the fields J^0_t, J^1_t and the Jacobians |\det D\varphi_{t,1}|; the kernel applied by FFT; a disk-to-C test pair and two smoothed synthetic “brain slices” built from ellipses.
   • Success criterion: relative mismatch below 10 % on the disk-to-C pair with a positive Jacobian everywhere; agreement between the explicit gradient and an automatic-differentiation gradient (Code 6.6) to 1 % at a random point.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 3 and 6.
2. [Solo · ≈ 20 h] Paper study: varifolds for unoriented shapes. Reproduce the main claim of Charon and Trouvé (2013): varifold matching does not need consistent orientations, while current matching does.
   • Work: extend Code 5.7 with a current data term (Proposition 5.4); flip the orientation of a random fraction of the segments of the target; repeat on a triangulated surface in 3D (a sphere matched to an ellipsoid with a bump).
   • Success criterion: a curve of final varifold distance against the fraction of flipped elements, from 0 to 50 %, flat for varifolds and degrading for currents.
   • Deliverables: code, report, talk. Chapters 4 and 5.
3. [Solo · ≈ 20 h] Paper study: statistics via initial momenta. Reproduce the approach of Vaillant, Miller, Younes and Trouvé (2004) on synthetic data with a known answer.
   • Work: generate populations as in Code 8.4 with two or three known modes; match the true template to each subject (Code 4.7); compute the PCA of the momenta with the K(\bar q) inner product; compare with Kendall tangent PCA.
   • Success criterion: the angle between the estimated and the true two-dimensional mode subspaces below 10° at the noise level of Code 8.4; a curve of this angle against the noise.
   • Deliverables: code, report, talk. Chapters 1, 4 and 8.

Group projects

1. [Group of 2–3 · ≈ 30 h each] A surface atlas with varifolds. Build the template-estimation pipeline of Joshi, Davis, Jomier and Gerig (2004), with the varifold data term of Charon and Trouvé (2013), for synthetic 3D surfaces.
   • Data: 20 to 40 triangulated surfaces (10^3 to 5 \times 10^3 faces) obtained by shooting a bent ellipsoid with random momenta in two or three modes, then remeshing each one differently.
   • Split: one member writes the 3D shooting and the varifold term, one the alternating atlas estimation, one the statistics (PCA of momenta) and the evaluation.
   • Success criterion: the estimated template within a varifold distance of 5 % of the population spread from the generating template; the generating modes recovered as the first principal modes.
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 4, 5 and 8.
2. [Group of 2–3 · ≈ 30 h each] Registration of diffusion tensor fields. Compare the reorientation strategies of Alexander, Pierpaoli, Basser and Gee (2001) inside an LDDMM registration.
   • Data: synthetic 2D tensor fields with curved fibre bundles (tensors elongated along the tangent of a family of curves), deformed by known diffeomorphisms; noise added on the eigenvalues.
   • Split: one member builds the phantoms and the ground truth; one implements the registration of the tensor fields with a Frobenius and an affine-invariant data term (Proposition 2.5); one implements the full action, finite strain and preservation of principal direction, and evaluates them.
   • Success criterion: a table of the mean angle between registered and true fibre directions, for the three strategies and two data terms. The best combination must bring this angle below half of its value before registration, and the ranking of the strategies must be explained on the data with Example 6.2.
   • Deliverables: code, report, talk. Chapters 2 and 6.
3. [Group of 2–3 · ≈ 30 h each] Lesions and metamorphosis. Compare LDDMM and image metamorphosis (Trouvé and Younes, 2005) for matching a healthy template to images with a lesion.
   • Data: synthetic 2D “brain” phantoms made of nested smoothed ellipses, deformed, with a bright or dark lesion added in some of them.
   • Split: one member extends Code 6.6, one extends Code 7.6 to 2D, one builds the phantoms, chooses \sigma and runs the evaluation.
   • Success criterion: a curve of the error on the recovered deformation (known for the phantoms), measured outside the lesion, against the lesion size, for LDDMM and for metamorphosis. For lesions up to 10 % of the phantom area, metamorphosis should keep this error within 20 % of the lesion-free case.
   • Deliverables: code, report, talk. Chapters 6 and 7.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet allowed, no electronic devices, four exercises.

1. Question (4 points). Explain why large deformation models use the flow of a time-dependent velocity field rather than the map x \mapsto x + v(x). State the regularity of the flow and the formula for its Jacobian determinant.

   Expected: x + v(x) is a diffeomorphism when \|Dv\|_\infty < 1, but folds where \det(I + Dv) < 0; the flow of an L^2([0, 1], V) velocity with V admissible is a C^1 diffeomorphism (Theorem 3.3); the Jacobian is the exponential of the integral of the divergence along the trajectory, hence positive (Proposition 0.7).

2. Question (6 points). Gaussian kernel of width \sigma, two landmarks at distance \sigma/2. (a) Compute the cost of moving one landmark with velocity u while the other stays still. (b) Same question for landmark metamorphosis with \sigma_{\text{meta}}^2 = 1. (c) Compute the affine-invariant distance between the identity and \operatorname{diag}(4, 1), and the determinants of their Euclidean and affine-invariant midpoints.

   Expected: (a) \kappa = e^{-1/4} = 0.779 and the cost is |u|^2 / (1 - \kappa^2) = 2.541\, |u|^2. (b) 2 / (4 - \kappa^2) = 0.589 times |u|^2. (c) The distance is \log 4 = 1.386; the determinants are 2.5 and 2.

3. Question (5 points). Let a finite-dimensional group act on Q with infinitesimal action A_q onto T_q Q. Prove that the projected metric is |\eta|_q^2 = \eta^\top K(q)^{-1} \eta with K(q) = A_q A_q^*. Deduce the form of the optimal velocity field in landmark matching.

   Expected: the minimal-norm solution of A_q \xi = \eta lies in the orthogonal of \ker A_q, which is the range of A_q^*. Writing \xi = A_q^* p gives K p = \eta and |\xi|^2 = p^\top K p.

   For landmarks, A_q^* p = \sum_j K(\cdot, q_j) p_j: the optimal velocity is a combination of kernels at the landmarks (Theorems 2.8 and 4.2).

4. Question (5 points). Write the LDDMM energy for images and its gradient. Explain why optimal momenta are scalar fields times the image gradient, and what goes wrong when the target contains a structure absent from the template; describe how metamorphosis fixes it.

   Expected: Definition 6.4 and Theorem 6.5; Proposition 6.3 (velocities along level lines are invisible, so the optimal momentum lies in the range of A_I^*); a new structure is out of the orbit, so the distance is infinite and the optimizer distorts the anatomy; metamorphosis adds the intensity residual Z with weight 1/\sigma^2, giving a finite distance (Proposition 7.2).

Grading rubric (proposal)

References

Official references

Neither the master-mva.com course page nor the Université Paris-Saclay course detail lists references for this course. The textbook of Younes (2010) and the papers of the instructors below cover most of the official syllabus.

Textbooks

• L. Younes, Shapes and Diffeomorphisms, Applied Mathematical Sciences 171, Springer, 2010 (2nd ed. 2019).
• U. Grenander and M. I. Miller, Pattern Theory: From Representation to Inference, Oxford University Press, 2007.
• I. L. Dryden and K. V. Mardia, Statistical Shape Analysis, with Applications in R, 2nd ed., Wiley, 2016.
• D. G. Kendall, D. Barden, T. K. Carne and H. Le, Shape and Shape Theory, Wiley, 1999.
• A. Srivastava and E. Klassen, Functional and Shape Data Analysis, Springer, 2016.
• X. Pennec, S. Sommer and T. Fletcher (eds.), Riemannian Geometric Statistics in Medical Image Analysis, Academic Press, 2020.
• M. P. do Carmo, Riemannian Geometry, Birkhäuser, 1992.
• J. M. Lee, Introduction to Riemannian Manifolds, 2nd ed., Graduate Texts in Mathematics 176, Springer, 2018.
• D. D. Holm, T. Schmah and C. Stoica, Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions, Oxford University Press, 2009.
• R. Bhatia, Positive Definite Matrices, Princeton University Press, 2007.

Foundations: kernels, Riemannian geometry, groups of diffeomorphisms (ch. 0–3)

• N. Aronszajn, “Theory of reproducing kernels”, Transactions of the American Mathematical Society 68(3), 1950.
• C. J. Atkin, “The Hopf–Rinow theorem is false in infinite dimensions”, Bulletin of the London Mathematical Society 7(3), 1975.
• D. G. Kendall, “Shape manifolds, Procrustean metrics, and complex projective spaces”, Bulletin of the London Mathematical Society 16(2), 1984.
• J. C. Gower, “Generalized Procrustes analysis”, Psychometrika 40(1), 1975.
• X. Pennec, P. Fillard and N. Ayache, “A Riemannian framework for tensor computing”, International Journal of Computer Vision 66(1), 2006.
• V. Arsigny, P. Fillard, X. Pennec and N. Ayache, “Geometric means in a novel vector space structure on symmetric positive-definite matrices”, SIAM Journal on Matrix Analysis and Applications 29(1), 2007.
• V. I. Arnold, “Sur la géométrie différentielle des groupes de Lie de dimension infinie et ses applications à l’hydrodynamique des fluides parfaits”, Annales de l’Institut Fourier 16(1), 1966.
• D. D. Holm, J. E. Marsden and T. S. Ratiu, “The Euler–Poincaré equations and semidirect products with applications to continuum theories”, Advances in Mathematics 137(1), 1998.
• A. Trouvé, “Diffeomorphisms groups and pattern matching in image analysis”, International Journal of Computer Vision 28(3), 1998.
• P. Dupuis, U. Grenander and M. I. Miller, “Variational problems on flows of diffeomorphisms for image matching”, Quarterly of Applied Mathematics 56(3), 1998.
• M. I. Miller, A. Trouvé and L. Younes, “On the metrics and Euler–Lagrange equations of computational anatomy”, Annual Review of Biomedical Engineering 4, 2002.
• M. I. Miller, A. Trouvé and L. Younes, “Geodesic shooting for computational anatomy”, Journal of Mathematical Imaging and Vision 24(2), 2006.

Landmarks, point sets, curves and surfaces (ch. 4–5)

• S. C. Joshi and M. I. Miller, “Landmark matching via large deformation diffeomorphisms”, IEEE Transactions on Image Processing 9(8), 2000.
• M. Micheli, P. W. Michor and D. Mumford, “Sectional curvature in terms of the cometric, with applications to the Riemannian manifolds of landmarks”, SIAM Journal on Imaging Sciences 5(1), 2012.
• L. Risser, F.-X. Vialard, R. Wolz, M. Murgasova, D. D. Holm and D. Rueckert, “Simultaneous multi-scale registration using large deformation diffeomorphic metric mapping”, IEEE Transactions on Medical Imaging 30(10), 2011.
• J. Glaunès, A. Trouvé and L. Younes, “Diffeomorphic matching of distributions: a new approach for unlabelled point-sets and sub-manifolds matching”, IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2004.
• M. Vaillant and J. Glaunès, “Surface matching via currents”, Information Processing in Medical Imaging (IPMI), LNCS 3565, 2005.
• J. Glaunès, A. Qiu, M. I. Miller and L. Younes, “Large deformation diffeomorphic metric curve mapping”, International Journal of Computer Vision 80(3), 2008.
• N. Charon and A. Trouvé, “The varifold representation of nonoriented shapes for diffeomorphic registration”, SIAM Journal on Imaging Sciences 6(4), 2013.
• P. Roussillon and J. Glaunès, “Kernel metrics on normal cycles and application to curve matching”, SIAM Journal on Imaging Sciences 9(4), 2016.
• I. Kaltenmark, B. Charlier and N. Charon, “A general framework for curve and surface comparison and registration with oriented varifolds”, IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2017.

Images, tensor fields and metamorphoses (ch. 6–7)

• M. F. Beg, M. I. Miller, A. Trouvé and L. Younes, “Computing large deformation metric mappings via geodesic flows of diffeomorphisms”, International Journal of Computer Vision 61(2), 2005.
• D. C. Alexander, C. Pierpaoli, P. J. Basser and J. C. Gee, “Spatial transformations of diffusion tensor magnetic resonance images”, IEEE Transactions on Medical Imaging 20(11), 2001.
• Y. Cao, M. I. Miller, R. L. Winslow and L. Younes, “Large deformation diffeomorphic metric mapping of vector fields”, IEEE Transactions on Medical Imaging 24(9), 2005.
• M. I. Miller and L. Younes, “Group actions, homeomorphisms, and matching: a general framework”, International Journal of Computer Vision 41(1–2), 2001.
• A. Trouvé and L. Younes, “Metamorphoses through Lie group action”, Foundations of Computational Mathematics 5(2), 2005.
• D. D. Holm, A. Trouvé and L. Younes, “The Euler–Poincaré theory of metamorphosis”, Quarterly of Applied Mathematics 67(4), 2009.

Statistics and computational anatomy (ch. 8)

• H. Karcher, “Riemannian center of mass and mollifier smoothing”, Communications on Pure and Applied Mathematics 30(5), 1977.
• B. Afsari, “Riemannian L^p center of mass: existence, uniqueness, and convexity”, Proceedings of the American Mathematical Society 139(2), 2011.
• P. T. Fletcher, C. Lu, S. M. Pizer and S. Joshi, “Principal geodesic analysis for the study of nonlinear statistics of shape”, IEEE Transactions on Medical Imaging 23(8), 2004.
• M. Vaillant, M. I. Miller, L. Younes and A. Trouvé, “Statistics on diffeomorphisms via tangent space representations”, NeuroImage 23(Suppl. 1), 2004.
• S. Joshi, B. Davis, M. Jomier and G. Gerig, “Unbiased diffeomorphic atlas construction for computational anatomy”, NeuroImage 23(Suppl. 1), 2004.
• S. Allassonnière, Y. Amit and A. Trouvé, “Towards a coherent statistical framework for dense deformable template estimation”, Journal of the Royal Statistical Society: Series B 69(1), 2007.
• S. Allassonnière, E. Kuhn and A. Trouvé, “Construction of Bayesian deformable models via a stochastic approximation algorithm: a convergence study”, Bernoulli 16(3), 2010.
• J. G. Csernansky, S. Joshi, L. Wang, J. W. Haller, M. Gado, J. P. Miller, U. Grenander and M. I. Miller, “Hippocampal morphometry in schizophrenia by high dimensional brain mapping”, Proceedings of the National Academy of Sciences 95(19), 1998.
• S. Durrleman, X. Pennec, A. Trouvé and N. Ayache, “Statistical models of sets of curves and surfaces based on currents”, Medical Image Analysis 13(5), 2009.
• S. Durrleman, X. Pennec, A. Trouvé, J. Braga, G. Gerig and N. Ayache, “Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data”, International Journal of Computer Vision 103(1), 2013.
• M. Niethammer, Y. Huang and F.-X. Vialard, “Geodesic regression for image time-series”, Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2011.
• P. T. Fletcher, “Geodesic regression and the theory of least squares on Riemannian manifolds”, International Journal of Computer Vision 105(2), 2013.

Research directions (ch. 9)

• P. W. Michor and D. Mumford, “Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms”, Documenta Mathematica 10, 2005.
• P. W. Michor and D. Mumford, “Riemannian geometries on spaces of plane curves”, Journal of the European Mathematical Society 8(1), 2006.
• L. Younes, “Computable elastic distances between shapes”, SIAM Journal on Applied Mathematics 58(2), 1998.
• A. Srivastava, E. Klassen, S. H. Joshi and I. H. Jermyn, “Shape analysis of elastic curves in Euclidean spaces”, IEEE Transactions on Pattern Analysis and Machine Intelligence 33(7), 2011.
• M. Bauer, M. Bruveris and P. W. Michor, “Overview of the geometries of shape spaces and diffeomorphism groups”, Journal of Mathematical Imaging and Vision 50(1–2), 2014.
• L. Devilliers, S. Allassonnière, A. Trouvé and X. Pennec, “Template estimation in computational anatomy: Fréchet means in top and quotient spaces are not consistent”, SIAM Journal on Imaging Sciences 10(3), 2017.
• B. Gris, S. Durrleman and A. Trouvé, “A sub-Riemannian modular framework for diffeomorphism-based analysis of shape ensembles”, SIAM Journal on Imaging Sciences 11(1), 2018.
• S. Arguillère, E. Trélat, A. Trouvé and L. Younes, “Shape deformation analysis from the optimal control viewpoint”, Journal de Mathématiques Pures et Appliquées 104(1), 2015.
• A. Trouvé and F.-X. Vialard, “Shape splines and stochastic shape evolutions: a second order point of view”, Quarterly of Applied Mathematics 70(2), 2012.
• A. Arnaudon, D. D. Holm and S. Sommer, “A geometric framework for stochastic shape analysis”, Foundations of Computational Mathematics 19(3), 2019.
• J. Feydy, B. Charlier, F.-X. Vialard and G. Peyré, “Optimal transport for diffeomorphic registration”, Medical Image Computing and Computer-Assisted Intervention (MICCAI), 2017.
• A. Genevay, G. Peyré and M. Cuturi, “Learning generative models with Sinkhorn divergences”, International Conference on Artificial Intelligence and Statistics (AISTATS), 2018.
• J. Feydy, T. Séjourné, F.-X. Vialard, S. Amari, A. Trouvé and G. Peyré, “Interpolating between optimal transport and MMD using Sinkhorn divergences”, International Conference on Artificial Intelligence and Statistics (AISTATS), 2019.
• B. Charlier, J. Feydy, J. A. Glaunès, F.-D. Collin and G. Durif, “Kernel operations on the GPU, with autodiff, without memory overflows”, Journal of Machine Learning Research 22(74), 2021.
• A. Bône, M. Louis, B. Martin and S. Durrleman, “Deformetrica 4: an open-source software for statistical shape analysis”, Shape in Medical Imaging (ShapeMI), MICCAI workshop, LNCS 11167, 2018.
• N. Miolane, N. Guigui, A. Le Brigant and coauthors, “Geomstats: a Python package for Riemannian geometry in machine learning”, Journal of Machine Learning Research 21(223), 2020.
• X. Yang, R. Kwitt, M. Styner and M. Niethammer, “Quicksilver: fast predictive image registration – a deep learning approach”, NeuroImage 158, 2017.
• G. Balakrishnan, A. Zhao, M. R. Sabuncu, J. Guttag and A. V. Dalca, “VoxelMorph: a learning framework for deformable medical image registration”, IEEE Transactions on Medical Imaging 38(8), 2019.
• A. V. Dalca, G. Balakrishnan, J. Guttag and M. R. Sabuncu, “Unsupervised learning of probabilistic diffeomorphic registration for images and surfaces”, Medical Image Analysis 57, 2019.