Deformable Models and Geodesic Methods for Image Analysis

Contents

• Introduction and overview
• 0. Review: calculus of variations, plane curves and discretisation
• 1. Active contours: snakes and balloons
• 2. Finite differences, finite elements and deformable surfaces
• 3. Level sets and geodesic active contours
• 4. Active regions and shape priors
• 5. Minimal paths and geodesics
• 6. The eikonal equation and Fast Marching
• 7. Metrics: 2D, 3D, surfaces, anisotropic, space and radius
• 8. Geodesic remeshing and applications
• 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 who knows multivariate calculus and some Python, and who wants to extract contours, surfaces and thin structures from medical or natural images with methods whose behaviour can be proved.

Empirical first: put a closed elastic curve of radius 8 pixels inside a noisy disk of radius 30 and let it minimise its energy. It never reaches the edge: it only shrinks a little, to a mean radius of 7.2, because the edge is too far away to be felt (Code 1.11). Push it from inside with a small pressure and it stops at a radius of 30.45; push a little harder and it bursts through the edge and leaves the image.

Shortest paths have a similar surprise: on a 4-connected pixel grid, the shortest path between two pixels overestimates the diagonal distance by 41 % (Code 5.11). This course explains both phenomena and builds the methods that cure them.

The course has two halves, as the official organisation announces. The first half treats deformable models: curves and surfaces that move to minimise an energy combining regularity and image fit, from parametric snakes and balloons (ch. 1–2) to level sets, geodesic active contours, active regions and shape priors (ch. 3–4).

The second half treats geodesic methods: the globally minimal path between two points for a metric built from the image, computed by solving the eikonal equation with the Fast Marching algorithm (ch. 5–6). Metrics on volumes, surfaces, anisotropic and space–radius domains (ch. 7) lead to geodesic remeshing and to the medical applications of the course (ch. 8).

Image models and denoising are developed in The Machine Intelligence of Images and Image Denoising. Shapes as points of a space are the subject of Geometry and Shape Spaces. Medical image analysis is also treated in Medical Image Analysis (semester 1) and Deep Learning for Medical Imaging (semester 2).

What you will be able to do

• Derive the Euler–Lagrange equation of a curve or surface energy and turn it into a gradient flow.
• Implement a snake with a stable semi-implicit scheme, explain its spectrum, and add a balloon force or a distance-based potential to enlarge its capture range.
• Compare finite-difference and finite-element discretisations of an active contour, and write the energy of a deformable surface.
• Represent an evolving front by a level set, derive the geodesic active contour, and discretise it with upwind schemes that allow topology changes.
• Segment images with region-based energies (Chan–Vese, region competition) and add a statistical shape prior.
• Define a geodesic distance from an image, prove that it solves the eikonal equation, and extract a minimal path by back-propagation.
• Implement Fast Marching, prove its causality and use competing fronts to build geodesic Voronoi diagrams and to link pairs of points.
• Design isotropic, anisotropic and space–radius metrics for edges, lines and vessels, in 2D, in 3D and on surfaces.
• Sample and remesh a domain or a surface by farthest-point geodesic sampling, and extract centrelines, trees and tubes from medical images.
• Read a research paper on deformable models or minimal paths, reimplement it and evaluate it, as the course project asks.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com sets the objective of giving a large overview of the methods and algorithms of deformable models for image segmentation, illustrated by concrete applications. It links to a presentation of the course, which these notes do not reproduce.

The same page lists eleven themes (“Thèmes abordés”) followed by the heading of the instructors list, and says that the course is “about half on deformable models and half on geodesic methods”. The Université Paris-Saclay course detail gives the volume (21 h of lectures, 9 h of labs) and the nature of the assessment, but no programme items. Each official item is matched below.

Beyond the syllabus, this page adds a review chapter (ch. 0), a research chapter for self-study on curvature-penalised minimal paths, convex relaxations and learned deformable models (ch. 9), and a written-exam proposal.

Where this course fits in the programme

Official prerequisite: “Admission to MVA sufficient”. The course page on master-mva.com lists Laurent Cohen (CNRS and Université Paris-Dauphine PSL) as the instructor. The course carries the Santé (health) track label.

• Refreshers — Linear Algebra, Analysis and Optimization for gradients, integration by parts and the spectral theorem; Scientific Python, NumPy and PyTorch for the labs; Probability and Statistics for the region models of ch. 4.
• Semester 1 courses that prepare it — Convex Optimization and Applications in Machine Learning (gradient flows, convex relaxations), The Machine Intelligence of Images (image models), Sub-Pixel Image Processing (interpolation), Medical Image Analysis and Geometric Data Analysis.
• Related courses of semester 2 — Geometry and Shape Spaces (shape metrics, deformations), Inverse Problems and Imaging (variational models), Point Clouds and 3D Modelling (surfaces), Deep Learning for Medical Imaging and Representation Learning for Computer Vision and Medical Imaging.
• Afterwards — the Research Internship, for which the paper-study project of this course is a rehearsal.

Course outline

• Review: calculus of variations, plane curves and discretisation — Euler–Lagrange equations, gradient flows, curvature, stability of finite differences, Dijkstra’s algorithm.
• Active contours: snakes and balloons — the snake energy, its semi-implicit evolution, and the forces that enlarge its capture range.
• Finite differences, finite elements and deformable surfaces — two discretisations compared, and the extension from curves to surfaces in 3D images.
• Level sets and geodesic active contours — implicit fronts, topology changes, weighted length, upwind schemes and competing fronts.
• Active regions and shape priors — region-based energies, region competition, and statistical models of shape.
• Minimal paths and geodesics — globally optimal curves between two points, the surface of minimal action and back-propagation.
• The eikonal equation and Fast Marching — viscosity solutions, the upwind scheme, the Fast Marching algorithm and front competition.
• Metrics: 2D, 3D, surfaces, anisotropic, space and radius — how to design the metric for edges, lines, vessels and surfaces.
• Geodesic remeshing and applications — farthest-point sampling, geodesic meshes, vessel trees, virtual endoscopy and surface segmentation.
• Research directions and open problems — curvature-penalised paths, global minimisers, heat-based distances and learned deformable models (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. The organisation announced on master-mva.com differs: 11 sessions of 3 hours, three of which include 2 hours of lab, which would make 33 h. The plan below keeps the 30 h of the programme table, with seven lecture sessions and three lab sessions of 3 h; the split into sessions is a proposal.

Level 1 = foundations, 2 = core models and algorithms, 3 = advanced metrics and applications close to research; the level never decreases from one session to the next. Sessions 1–5 cover deformable models, sessions 6–10 geodesic methods.

Self-study, outside the 30 h (about 4 h in total): ch. 9 Research directions and open problems — curvature-penalised and Finsler minimal paths, global minimisers, heat-based distances, learned deformable models (≈ 4 h, with Ex. 18–19).

How to study this module

• Read ch. 0 first even if the material looks familiar: the Euler–Lagrange equation, the first variation of length and the stability condition of finite differences return in every chapter.
• Keep the two halves connected: the geodesic active contour of ch. 3 and the minimal path of ch. 5 minimise the same weighted length, one locally by a flow, the other globally by the eikonal equation.
• After each session, run the code of the chapter and change one parameter (the pressure, the edge-stopping constant, the metric) to see what breaks.
• Do the exercises listed in the session plan before the next session; the ★★★ ones are good starting points for the project.
• Python environment: NumPy and SciPy (scipy.ndimage for filtering and distance transforms, scipy.sparse.csgraph for Dijkstra); every code block runs in a few seconds on a laptop with synthetic images. Readers who prefer Matlab can translate the array code directly.
• Start choosing a paper for the project around session 6: the list of solo subjects gives the scale to aim at.

0. Review: calculus of variations, plane curves and discretisation

Every model of this course is an energy on curves or surfaces, minimised by a flow that a computer solves on a grid. This review states the tools used throughout: the Euler–Lagrange equation, gradient flows, the geometry of plane curves, the stability of finite differences, and Dijkstra’s algorithm.

Energies and their first variation

We work with closed curves parametrised on [0, 1], so every function below is 1-periodic and integrations by parts have no boundary terms.

Plane curves

Finite differences and stability

On a grid of step h, derivatives become differences. For a periodic sequence u_i \approx u(ih): D^+ u_i = \frac{u_{i+1} - u_i}{h}, \qquad D^- u_i = \frac{u_i - u_{i-1}}{h}, \qquad D^2 u_i = \frac{u_{i+1} - 2 u_i + u_{i-1}}{h^2} .

import numpy as np

# Heat equation x_t = x_ss on a closed curve (periodic), n nodes, step h = 1/n
n = 64
h = 1.0 / n
s = np.arange(n) * h
rng = np.random.default_rng(0)
x0 = np.cos(2 * np.pi * s) + 0.05 * rng.standard_normal(n)   # noisy circle coordinate

D2 = (np.roll(np.eye(n), 1, 0) - 2 * np.eye(n) + np.roll(np.eye(n), -1, 0)) / h**2

def explicit(x, tau, steps):
    for _ in range(steps):
        x = x + tau * D2 @ x
    return x

def implicit(x, tau, steps):
    M = np.linalg.inv(np.eye(n) - tau * D2)
    for _ in range(steps):
        x = M @ x
    return x

for r in [0.45, 0.55]:
    tau = r * h**2
    print(f"explicit, tau = {r} h^2: max |x| after 400 steps = {np.abs(explicit(x0, tau, 400)).max():.3g}")
print(f"implicit, tau = 50 h^2: max |x| after 400 steps = {np.abs(implicit(x0, 50 * h**2, 400)).max():.3g}")
lam = np.linalg.eigvalsh(-D2)
print(f"largest eigenvalue of -D2: {lam.max():.1f} = 4/h^2 = {4 / h**2:.1f}")

Images, edges and shortest paths

1. Active contours: snakes and balloons

How can a computer draw the outline of an organ in a noisy image? The active contour answer is to place an elastic curve near the object and let it slide down an energy that is low on edges. This chapter defines that energy, computes its evolution stably, and fixes its main weakness, the small capture range.

The snake energy

Discretisation and stability

The capture range and the balloon force

In a discrete implementation, the normalisation by the local norm makes the noise of flat regions as strong as a real edge. Code 1.11 therefore divides \nabla P by its maximum over the image, so that the image force equals k only on the strongest edges and fades in flat regions.

import numpy as np
from scipy import ndimage as ndi

# Synthetic image: a bright disk of radius 30 on a dark background, with noise
N, R_true = 128, 30.0
yy, xx = np.mgrid[0:N, 0:N].astype(float)
c = np.array([64.0, 64.0])
rng = np.random.default_rng(1)
I = (np.hypot(xx - c[0], yy - c[1]) < R_true).astype(float) + 0.1 * rng.standard_normal((N, N))

# Potential P = -|grad(G_sigma * I)|^2 and its gradient
Is = ndi.gaussian_filter(I, 2.0)
P = -(ndi.sobel(Is, 1) ** 2 + ndi.sobel(Is, 0) ** 2)
Py, Px = np.gradient(P)

def sample(F, v):                       # bilinear interpolation at points v = (x, y)
    return ndi.map_coordinates(F, [v[:, 1], v[:, 0]], order=1)

# Closed snake: n nodes, internal matrix A (pentadiagonal, circulant)
n, a, b, tau = 80, 0.05, 0.01, 1.0
row = np.zeros(n); row[[0, 1, 2, -2, -1]] = [2 * a + 6 * b, -a - 4 * b, b, b, -a - 4 * b]
A = np.array([np.roll(row, k) for k in range(n)])
M = np.linalg.inv(np.eye(n) + tau * A)  # semi-implicit step, computed once

theta = 2 * np.pi * np.arange(n) / n
v = c + 8.0 * np.c_[np.cos(theta), np.sin(theta)]   # start inside, radius 8

def normals(v):                          # outward unit normals of a counter-clockwise curve
    t = np.roll(v, -1, 0) - np.roll(v, 1, 0)
    nrm = np.c_[t[:, 1], -t[:, 0]]
    return nrm / np.linalg.norm(nrm, axis=1, keepdims=True)

def run(k1, k, steps=300):
    v_ = v.copy()
    for _ in range(steps):
        g = np.c_[sample(Px, v_), sample(Py, v_)]
        gn = np.sqrt(Px**2 + Py**2).max()          # global normalisation of the image force
        F = k1 * normals(v_) - k * g / gn  # balloon force + image force
        v_ = M @ (v_ + tau * F)
    return np.linalg.norm(v_ - c, axis=1)

for k1 in [0.0, 0.5, 1.2]:
    r = run(k1=k1, k=1.0)
    print(f"k1 = {k1}: radius mean {r.mean():6.2f}, min {r.min():6.2f}, max {r.max():6.2f}  (true {R_true})")
lam = np.linalg.eigvalsh(A)
print(f"eigenvalues of A in [{lam.min():.2e}, {lam.max():.3f}], 4a + 16b = {4 * a + 16 * b:.3f}")

2. Finite differences, finite elements and deformable surfaces

Finite differences sample the snake at its nodes; finite elements represent it as a continuous curve built from local basis functions. This chapter compares the two on the same equation, then extends active contours from curves in 2D images to surfaces in 3D images.

The weak form and finite elements

Multiply the snake flow of Proposition 1.2 by a test function w and integrate by parts on the closed curve. The equation holds in the weak form \int_0^1 \big( \langle v_t, w \rangle + \alpha \langle v', w' \rangle + \beta \langle v'', w'' \rangle \big) \, ds = \int_0^1 \langle F(v), w \rangle \, ds for every periodic test function w, where F = -\nabla P. It needs one derivative less than the strong form.

import numpy as np
from scipy.linalg import eigh

# Membrane term -v'' on a closed curve of length 1, n nodes, h = 1/n.
n = 32
h = 1.0 / n
I = np.eye(n)
S = np.roll(I, 1, 1) + np.roll(I, -1, 1)              # neighbours i-1 and i+1 (periodic)
K = (2 * I - S) / h                                   # P1 stiffness:  int phi_i' phi_j'
M = h * (4 * I + S) / 6                               # P1 consistent mass: int phi_i phi_j
M_lumped = h * I                                      # lumped mass = finite differences

lam_fe = eigh(K, M, eigvals_only=True)
lam_fd = eigh(K, M_lumped, eigvals_only=True)
for k in [1, 4, 8]:
    exact = (2 * np.pi * k) ** 2
    # each nonzero frequency appears twice (cos and sin); index 2k-1 picks frequency k
    print(f"k = {k}: exact {exact:8.1f} | FE {lam_fe[2 * k - 1]:8.1f} | FD {lam_fd[2 * k - 1]:8.1f}")

What finite elements bring to active contours is less this accuracy than the representation. The curve is defined between the nodes, and the image force enters through the integrals \int F(v) \varphi_i along each element, so an edge crossing an element is felt even between two nodes.

This allows fewer nodes, and it extends directly to surfaces meshed by triangles. Cohen and Cohen (1993) develop finite-element active contours and balloons for 2D and 3D images, an official reference of the course.

Deformable surfaces

3. Level sets and geodesic active contours

A parametric snake cannot split around two objects or merge with another snake. This chapter represents the moving curve as the zero level of a function on the whole image, which handles topology changes for free, and derives the geodesic active contour, the bridge between deformable models and the geodesic methods of the second half.

Implicit fronts

From geometric to geodesic active contours

The new term \langle \nabla g, \nabla \phi \rangle is an advection along -\nabla g, towards the bottom of the valley of g. It can hold the front on an edge where g is not zero, and pull it back after an overshoot, which the geometric model cannot do.

Numerical schemes

Two devices keep the cost and the accuracy under control. Re-initialisation periodically replaces \phi by the signed distance to its zero level, which keeps |\nabla \phi| \approx 1; computing a distance is an eikonal problem, solved by the methods of ch. 6. The narrow band method of Adalsteinsson and Sethian (1995) updates \phi only near the front.

import numpy as np
from scipy import ndimage as ndi

N = 128
yy, xx = np.mgrid[0:N, 0:N].astype(float)
truth = (np.hypot(xx - 40, yy - 64) < 15) | (np.hypot(xx - 88, yy - 64) < 15)
rng = np.random.default_rng(2)
I = truth + 0.1 * rng.standard_normal((N, N))

Is = ndi.gaussian_filter(I, 1.5)
grad2 = ndi.sobel(Is, 0) ** 2 + ndi.sobel(Is, 1) ** 2
g = 1.0 / (1.0 + grad2 / 0.3)                        # edge-stopping function, about 1 in flat areas
gy, gx = np.gradient(g)

def sdf(inside):                                     # signed distance, negative inside
    return ndi.distance_transform_edt(~inside) - ndi.distance_transform_edt(inside)

phi = sdf(np.hypot(xx - 64, yy - 64) < 55)           # one curve around both objects
nu, dt = 1.0, 0.2
for it in range(800):
    dxm = phi - np.roll(phi, 1, 1); dxp = np.roll(phi, -1, 1) - phi
    dym = phi - np.roll(phi, 1, 0); dyp = np.roll(phi, -1, 0) - phi
    # upwind |grad phi| for the shrinking speed (front moves inward: F = -nu g < 0)
    grad_minus = np.sqrt(np.minimum(dxm, 0)**2 + np.maximum(dxp, 0)**2
                         + np.minimum(dym, 0)**2 + np.maximum(dyp, 0)**2)
    # curvature term kappa |grad phi| with central differences
    py, px = np.gradient(phi)
    pyy, pyx = np.gradient(py); pxy, pxx = np.gradient(px)
    curv = (pxx * py**2 - 2 * px * py * pxy + pyy * px**2) / (px**2 + py**2 + 1e-8)
    # advection by grad g, upwinded
    adv = (np.maximum(gx, 0) * dxp + np.minimum(gx, 0) * dxm
           + np.maximum(gy, 0) * dyp + np.minimum(gy, 0) * dym)
    phi = phi + dt * (g * curv + nu * g * grad_minus + adv)
    if it % 20 == 19:
        phi = sdf(phi < 0)                           # re-initialisation
seg = phi < 0
labels, ncomp = ndi.label(seg)
iou = (seg & truth).sum() / (seg | truth).sum()
print(f"components: {ncomp}, IoU with the two disks: {iou:.3f}")

4. Active regions and shape priors

When an object differs from its background by its mean intensity or its texture, but not by a sharp edge, edge-based contours fail. This chapter moves the information from the boundary to the regions it encloses, then adds prior knowledge of the expected shape.

Region-based energies

To evolve the curve, write the energy with the level-set function \phi, negative inside, a smoothed Heaviside function H_\varepsilon and its derivative \delta_\varepsilon: H_\varepsilon(z) = \frac12 \Big( 1 + \frac{2}{\pi} \arctan \frac{z}{\varepsilon} \Big), \qquad \delta_\varepsilon(z) = \frac{\varepsilon}{\pi (\varepsilon^2 + z^2)} . The inside is H_\varepsilon(-\phi) \approx 1, and the length is \int \delta_\varepsilon(\phi) |\nabla \phi| \, dx.

Region competition and active regions

import numpy as np
from scipy import ndimage as ndi

N = 128
yy, xx = np.mgrid[0:N, 0:N].astype(float)
truth = ((np.abs(xx - 45) < 18) & (np.abs(yy - 60) < 25)) | (np.hypot(xx - 90, yy - 75) < 20)
rng = np.random.default_rng(3)
I = 0.3 + 0.4 * truth + 0.35 * rng.standard_normal((N, N))   # contrast 0.4, noise 0.35

def curvature(phi):
    py, px = np.gradient(phi)
    norm = np.sqrt(px**2 + py**2) + 1e-8
    return np.gradient(px / norm, axis=1) + np.gradient(py / norm, axis=0)

eps, mu, lam1, lam2, dt = 1.0, 0.5, 1.0, 1.0, 0.5
phi = np.hypot(xx - 64, yy - 64) - 30.0             # negative inside a central circle
for it in range(1000):
    inside = phi < 0
    c1, c2 = I[inside].mean(), I[~inside].mean()     # optimal constants (Proposition 4.3)
    delta = eps / (np.pi * (eps**2 + phi**2))        # smoothed Dirac
    force = mu * curvature(phi) + lam1 * (I - c1)**2 - lam2 * (I - c2)**2
    phi = phi + dt * delta * force / (np.abs(delta * force).max() + 1e-12)
    if it % 50 == 49:
        phi = ndi.distance_transform_edt(phi >= 0) - ndi.distance_transform_edt(phi < 0)

seg = phi < 0
iou = lambda s: (s & truth).sum() / (s | truth).sum()
print(f"c1 = {c1:.3f}, c2 = {c2:.3f} (true 0.7 and 0.3)")
print(f"Chan-Vese IoU: {iou(seg):.3f} | pixelwise threshold at 0.5: {iou(I > 0.5):.3f}")
sizes = np.bincount(ndi.label(seg)[0].ravel())[1:]
print(f"components: {sizes.size}, of which {np.sum(sizes > 50)} larger than 50 pixels")

Shape priors

When part of the object is occluded or invisible, as an organ next to a structure of the same intensity, image terms alone cannot recover it. A shape prior adds an energy that is low for the shapes seen in a training set.

import numpy as np
from scipy import ndimage as ndi

N = 96
yy, xx = np.mgrid[0:N, 0:N].astype(float)
def sdf(inside):                                       # signed distance, negative inside
    return ndi.distance_transform_edt(~inside) - ndi.distance_transform_edt(inside)
def ellipse(a, b):
    return ((xx - 48) / a) ** 2 + ((yy - 48) / b) ** 2 < 1

# Training set: 40 aligned ellipses with random semi-axes
rng = np.random.default_rng(6)
axes = rng.uniform(15, 35, size=(40, 2))
Phi = np.stack([sdf(ellipse(a, b)).ravel() for a, b in axes])
mean = Phi.mean(0)
U, s, Vt = np.linalg.svd(Phi - mean, full_matrices=False)
print("variance explained by the first 3 modes:", np.round(s[:3]**2 / (s**2).sum(), 3))

# A test shape with an occlusion: an ellipse with a bite taken out of it
true = ellipse(30, 20)
seen = true & ~(np.hypot(xx - 70, yy - 48) < 12)
w = Vt[:2] @ (sdf(seen).ravel() - mean)               # coordinates in the 2-mode shape space
proj = (mean + w @ Vt[:2]).reshape(N, N) < 0
iou = lambda s: (s & true).sum() / (s | true).sum()
print(f"IoU with the true ellipse: occluded {iou(seen):.3f}, projected on 2 modes {iou(proj):.3f}")

5. Minimal paths and geodesics

A snake finds the local minimum nearest to its initialisation. For an open curve with two given endpoints, such as a vessel between two points clicked by a radiologist, one can do much better: find the global minimum of a weighted length. This chapter, after Cohen and Kimmel (1997), defines the problem and shows how the solution is read off a single function, the surface of minimal action.

Geodesic distance

The surface of minimal action

Shortest paths on a pixel grid

import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import dijkstra

def grid_graph(P, offsets):
    """Graph on the pixels; edge u-v weighs |u - v| * (P(u) + P(v)) / 2."""
    H, W = P.shape
    idx = np.arange(H * W).reshape(H, W)
    rows, cols, w = [], [], []
    for dy, dx in offsets:
        a = idx[max(0, -dy):H - max(0, dy), max(0, -dx):W - max(0, dx)]
        b = idx[max(0, dy):H - max(0, -dy) or None, max(0, dx):W - max(0, -dx) or None]
        pa, pb = P.ravel()[a.ravel()], P.ravel()[b.ravel()]
        rows += [a.ravel()]; cols += [b.ravel()]; w += [np.hypot(dx, dy) * (pa + pb) / 2]
    r, c, w = map(np.concatenate, (rows, cols, w))
    return coo_matrix((w, (r, c)), shape=(H * W, H * W)).tocsr()

N4 = [(0, 1), (1, 0), (0, -1), (-1, 0)]
N8 = N4 + [(1, 1), (1, -1), (-1, 1), (-1, -1)]
N16 = N8 + [(1, 2), (2, 1), (-1, 2), (-2, 1), (1, -2), (2, -1), (-1, -2), (-2, -1)]

# (a) metrication error for the constant metric P = 1
n = 101
yy, xx = np.mgrid[0:n, 0:n]
euclid = np.hypot(xx - 50, yy - 50).ravel()
mask = euclid > 20
for name, off in [("4-neighbours", N4), ("8-neighbours", N8), ("16-neighbours", N16)]:
    d = dijkstra(grid_graph(np.ones((n, n)), off), indices=50 * n + 50)
    print(f"{name:14s}: max relative error {np.max(d[mask] / euclid[mask] - 1):.4f}")

# (b) a minimal path along a dark curve in a noisy image
H, W = 96, 128
yy, xx = np.mgrid[0:H, 0:W].astype(float)
curve = 48 + 25 * np.sin(2 * np.pi * xx[0] / W)
rng = np.random.default_rng(4)
I = np.clip(1.0 - (np.abs(yy - curve) < 1.5) + 0.3 * rng.standard_normal((H, W)), 0, None)
P = 0.05 + I                                        # cheap on the dark curve
G = grid_graph(P, N8)
src, dst = int(round(curve[0])) * W, int(round(curve[-1])) * W + W - 1
U, pred = dijkstra(G, indices=src, return_predecessors=True)
path, k = [], dst
while k != src:
    path.append(k); k = pred[k]
py, px = np.divmod(np.array(path), W)
print(f"path: {len(path)} pixels, mean |y - curve(x)| = {np.mean(np.abs(py - curve[px])):.2f} px")
print(f"action U(end) = {U[dst]:.1f}")

6. The eikonal equation and Fast Marching

Chapter 5 reduced minimal paths to one function, the minimal action U, which solves |\nabla U| = P. This chapter makes that statement precise, since U is not differentiable everywhere, and computes U on a grid in one pass with the Fast Marching algorithm.

Viscosity solutions

The upwind scheme and Fast Marching

In practice the error behaves like h away from singularities. Near a point source, where U is a cone, it decreases a little more slowly, as the next code shows.

import heapq
import numpy as np

def fast_marching(P, sources, h=1.0):
    """First-order Fast Marching for |grad U| = P on a grid of step h.
    Returns the action U and, for each point, the index of the source that reached it first."""
    H, W = P.shape
    U = np.full((H, W), np.inf)
    label = np.full((H, W), -1)
    alive = np.zeros((H, W), bool)
    heap = []
    for k, (i, j) in enumerate(sources):
        U[i, j], label[i, j] = 0.0, k
        heapq.heappush(heap, (0.0, i, j))
    while heap:
        u, i, j = heapq.heappop(heap)
        if alive[i, j]:
            continue
        alive[i, j] = True                                   # the smallest Trial value is final
        for a, b in ((i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)):
            if not (0 <= a < H and 0 <= b < W) or alive[a, b]:
                continue
            # smallest Alive neighbour along each axis (inf if none)
            ux = min(U[a, c] if 0 <= c < W and alive[a, c] else np.inf for c in (b - 1, b + 1))
            uy = min(U[c, b] if 0 <= c < H and alive[c, b] else np.inf for c in (a - 1, a + 1))
            f = h * P[a, b]
            if abs(ux - uy) >= f:                            # one-sided update
                new = min(ux, uy) + f
            else:                                            # two-sided (quadratic) update
                new = 0.5 * (ux + uy + np.sqrt(2 * f * f - (ux - uy) ** 2))
            if new < U[a, b]:
                U[a, b], label[a, b] = new, label[i, j]
                heapq.heappush(heap, (new, a, b))
    return U, label

# (a) first-order convergence for the Euclidean distance on [-1, 1]^2
for n in [51, 101, 201]:
    h = 2.0 / (n - 1)
    x = np.linspace(-1, 1, n)
    X, Y = np.meshgrid(x, x)
    U, _ = fast_marching(np.ones((n, n)), [(n // 2, n // 2)], h)
    print(f"n = {n:3d}, h = {h:.3f}: max error {np.max(np.abs(U - np.hypot(X, Y))):.4f}")

# (b) front competition between two sources in a non-constant metric
n = 121
x = np.linspace(-1, 1, n)
X, Y = np.meshgrid(x, x)
P = 1.0 + 2.0 * np.exp(-((X - 0.1) ** 2 + Y ** 2) / 0.05)   # a slow bump between the sources
p, q = (60, 20), (60, 100)
U, lab = fast_marching(P, [p, q], h=x[1] - x[0])
Up, _ = fast_marching(P, [p], h=x[1] - x[0])
border = (lab != np.roll(lab, 1, 1)) | (lab != np.roll(lab, 1, 0))
border[:, 0] = border[0, :] = False
print(f"d(p, q) = {Up[q]:.3f};  2 x min of U on the meeting line = {2 * U[border].min():.3f}")

Front competition and saddle points

7. Metrics: 2D, 3D, surfaces, anisotropic, space and radius

Fast Marching computes geodesics for any metric; what makes a minimal path follow a vessel rather than cut across the background is the metric. This chapter shows how to design it from the image, how to handle volumes and surfaces, and how two extensions, anisotropy and an extra radius dimension, capture the orientation and the width of tubular structures.

Isotropic metrics from images

Volumes and surfaces

Anisotropic metrics

Space and radius

import numpy as np
from scipy import ndimage as ndi
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import dijkstra

# A tube whose radius grows from 2 to 6 pixels along a sine-shaped centreline
H, W = 100, 128
yy, xx = np.mgrid[0:H, 0:W].astype(float)
xs = np.arange(10, 118)
cy = 50 + 15 * np.sin(2 * np.pi * xs / 120)
rad = 2 + 4 * (xs - 10) / 107
dist = np.full((H, W), np.inf)
for x0, y0, r0 in zip(xs, cy, rad):                     # union of disks along the centreline
    dist = np.minimum(dist, np.hypot(xx - x0, yy - y0) - r0)
rng = np.random.default_rng(5)
I = (dist < 0) + 0.25 * rng.standard_normal((H, W))

# Oriented flux on circles (2D version of Law and Chung, 2008):
# Q(x, r) = mean over the circle of radius r of grad(I_s)(x + r n) n^T, symmetrised.
# For a bright tube its smallest eigenvalue is most negative when r matches the radius.
Is = ndi.gaussian_filter(I, 1.0)
gy, gx = np.gradient(Is)
radii = np.arange(1.0, 8.01, 0.5)
ang = np.linspace(0, 2 * np.pi, 48, endpoint=False)
resp = []
for r in radii:
    Q = np.zeros((2, 2, H, W))
    for t in ang:
        n = np.array([np.cos(t), np.sin(t)])
        pts = [yy + r * n[1], xx + r * n[0]]
        g = [ndi.map_coordinates(gx, pts, order=1, mode="nearest"),
             ndi.map_coordinates(gy, pts, order=1, mode="nearest")]
        for i in range(2):
            for j in range(2):
                Q[i, j] += g[i] * n[j] / len(ang)
    a, b, c = Q[0, 0], (Q[0, 1] + Q[1, 0]) / 2, Q[1, 1]
    lam_min = (a + c) / 2 - np.sqrt(((a - c) / 2) ** 2 + b ** 2)
    resp.append(np.maximum(-lam_min, 0))
resp = np.stack(resp)                                            # shape (R, H, W)
P = np.exp(-8 * resp / resp.max()) + 0.01                       # metric on (r, y, x)

# Graph on (r, y, x): 8 spatial neighbours at fixed r, and r -> r +- 1
R = len(radii)
idx = np.arange(R * H * W).reshape(R, H, W)
rows, cols, wts = [], [], []
for dr, dy, dx in [(0, 0, 1), (0, 1, 0), (0, 1, 1), (0, 1, -1), (1, 0, 0)]:
    a = idx[:R - dr, max(0, -dy):H - dy, max(0, -dx):W - max(0, dx)]
    b = idx[dr:, dy:, max(0, dx):W - max(0, -dx) or None][:, :a.shape[1], :a.shape[2]]
    length = np.sqrt(dx**2 + dy**2 + (1.0 * dr)**2)              # a radius step of 0.5 costs like 1 pixel
    w = length * (P.ravel()[a.ravel()] + P.ravel()[b.ravel()]) / 2
    rows += [a.ravel(), b.ravel()]; cols += [b.ravel(), a.ravel()]; wts += [w, w]
G = coo_matrix((np.concatenate(wts), (np.concatenate(rows), np.concatenate(cols))),
               shape=(R * H * W,) * 2).tocsr()

start = idx[:, int(round(cy[0])), xs[0]]                         # all radii at the first end
U, pred, _ = dijkstra(G, indices=start, min_only=True, return_predecessors=True)
end = idx[:, int(round(cy[-1])), xs[-1]]
k = end[np.argmin(U[end])]
path = []
while k >= 0:
    path.append(k); k = pred[k]
pr, py, px = np.unravel_index(np.array(path), (R, H, W))
sel = (px >= xs[0]) & (px <= xs[-1])
err_c = np.abs(py[sel] - cy[px[sel] - xs[0]])
err_r = np.abs(radii[pr[sel]] - rad[px[sel] - xs[0]])
print(f"centreline error: mean {err_c.mean():.2f} px | radius error: mean {err_r.mean():.2f} px")
for x0 in [20, 60, 100]:
    j = np.flatnonzero(px == x0)[0]
    print(f"x = {x0:3d}: true radius {rad[x0 - 10]:.2f}, estimated {radii[pr[j]]:.1f}")

8. Geodesic remeshing and applications

Geodesic distances do more than link two points. Computed from many sources, they place samples evenly on a domain or a surface and build meshes adapted to a metric; computed from one root to many leaves, they extract trees. This chapter covers geodesic remeshing and the medical applications announced by the course.

Farthest-point sampling and geodesic meshes

import numpy as np
from scipy import ndimage as ndi
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import dijkstra

def graph(P, mask):
    """16-neighbour graph on the pixels of mask, edge weight |u - v| (P(u) + P(v)) / 2."""
    H, W = P.shape
    idx = np.arange(H * W).reshape(H, W)
    offs = [(0, 1), (1, 0), (1, 1), (1, -1), (1, 2), (2, 1), (1, -2), (2, -1)]
    r, c, w = [], [], []
    for dy, dx in offs:
        a = idx[:H - dy, max(0, -dx):W - max(0, dx)].ravel()
        b = a + dy * W + dx
        ok = mask.ravel()[a] & mask.ravel()[b]           # walls are >= 3 pixels thick: no jump
        a, b = a[ok], b[ok]
        wt = np.hypot(dx, dy) * (P.ravel()[a] + P.ravel()[b]) / 2
        r += [a, b]; c += [b, a]; w += [wt, wt]
    return coo_matrix((np.concatenate(w), (np.concatenate(r), np.concatenate(c))),
                      shape=(H * W, H * W)).tocsr()

# (a) farthest point sampling with the geodesic distance, in a room with a wall
H = W = 101
mask = np.ones((H, W), bool)
mask[:80, 49:52] = False                                 # a wall, open above row 80
G = graph(np.ones((H, W)), mask)
samples = [0]                                            # start at the corner (0, 0)
for k in range(1, 33):
    d = dijkstra(G, indices=samples, min_only=True)
    d[~mask.ravel()] = -1
    if k in (1, 2, 4, 8, 16, 32):
        print(f"k = {k:2d} samples: covering radius {d.max():6.2f},  radius x sqrt(k) = {d.max() * np.sqrt(k):6.1f}")
    samples.append(int(np.argmax(d)))

# (b) a centred path in a curved tube, as for virtual endoscopy
yy, xx = np.mgrid[0:H, 0:W]
tube = np.abs(np.hypot(xx - 50, yy - 100) - 60) < 8      # an arc of ring of width 16
d_wall = ndi.distance_transform_edt(tube)                # distance to the wall
ends = [np.flatnonzero((tube & (xx == x)).ravel()) for x in (2, 98)]
src_, dst = ends[0][len(ends[0]) // 2], ends[1][len(ends[1]) // 2]
for name, P in [("P = 1", np.ones((H, W))), ("P = 1/(d + 1)^2", 1.0 / (d_wall + 1.0) ** 2)]:
    U, pred = dijkstra(graph(P, tube), indices=src_, return_predecessors=True)
    k, path = dst, []
    while k != src_:
        path.append(k); k = pred[k]
    print(f"{name:16s}: smallest distance to the wall along the path = {d_wall.ravel()[path].min():.1f} px")

Trees, tubes, endoscopy and surfaces

9. Research directions and open problems

This chapter is for self-study. It points to four active directions that start from the methods of the course, with the papers where they are developed, and lists questions that remain open.

Curvature-penalised and Finsler minimal paths

A minimal path for an isotropic metric has no notion of direction: at a crossing of two vessels it may turn sharply onto the wrong branch, and between two structures it may take a shortcut through a corner. Penalising curvature fixes both, but the energy \int_C \big( 1 + \beta \, \kappa(s)^2 \big) \, P(C(s)) \, ds depends on the second derivative of the curve and is not a geodesic length in the image plane.

The way out is orientation lifting: a curve s \mapsto (x(s), \theta(s)) in the space of positions and orientations, constrained to move along its own direction. Curvature becomes the speed of \theta, and the energy becomes a length for a metric on the lifted space, sub-Riemannian or Finsler.

Bekkers, Duits, Mashtakov and Sanguinetti (2015) compute data-driven sub-Riemannian geodesics in this space with a PDE approach. Chen, Mirebeau and Cohen (2017) obtain global minimisers of a Finsler elastica model. Mirebeau (2018) designs Fast Marching schemes for several curvature-penalised models, among which the Reeds–Shepp car and the Euler–Mumford elastica.

Open problems: the lifted space has three dimensions in 2D and five in 3D, so memory and time are the bottleneck for volumes; and the choice of the curvature weight remains empirical.

Global minimisers for regions and surfaces

Chapter 4 minimised region energies by gradient flows, which find local minima. For two phases with fixed constants, a convex relaxation finds the global minimum.

Similar ideas lead to globally optimal minimal surfaces: discrete ones by graph cuts (Boykov and Kolmogorov, 2003), continuous ones by maximal flows (Appleton and Talbot, 2006). Pock, Cremers, Bischof and Chambolle (2009) extend convex relaxation to the full Mumford–Shah functional by lifting to a higher-dimensional space. Open problem: convex formulations with shape priors that remain tight.

Distances from heat

Varadhan (1967) proved that the heat kernel encodes the distance: for the heat kernel k_t of a Riemannian manifold, d(x, y) = \lim_{t \to 0} \sqrt{ -4 t \log k_t(x, y) } .

Crane, Weischedel and Wardetzky (2013) turn this into the heat method: solve the heat equation for a short time, normalise the gradient of the solution into a unit field X, then solve the Poisson equation \Delta \phi = \mathrm{div} X. Both steps are sparse linear systems whose factorisations can be reused for every new source, which makes the method fast on large meshes. Its accuracy near the cut locus and for strongly anisotropic metrics is still studied.

Deformable models and deep learning

Learning has entered deformable models in three ways.

• Learned energies. Marcos and coauthors (2018) predict the terms of a snake energy with a convolutional network and train the whole pipeline end to end.
• Learned contour evolution. Peng and coauthors (2020) deform an initial polygon with a network that applies circular convolutions to its vertices, for real-time instance segmentation.
• Geometric losses. Kervadec and coauthors (2019) train segmentation networks with a boundary loss computed from the distance map to the ground-truth boundary. Shit and coauthors (2021) propose clDice, a loss computed on skeletons that rewards the connectivity of tubular structures.

Open problems: a network does not guarantee that a vessel tree is connected or that an organ is simply connected, while minimal paths and level sets give such guarantees by construction. Learning the metric and keeping the geodesic extraction is one way to have both.

Questions to explore

• Crossings and bifurcations: how to extract a vascular tree where branches cross without merging, with orientation lifting or anisotropic metrics.
• Uncertainty: how confident is an extracted path, and how sensitive is it to the endpoints and the metric?
• Evaluation: overlap scores ignore topology; which measures reward correct connectivity and centrelines?
• Large volumes: partial propagation, multi-resolution and parallel solvers for the eikonal equation in 10^9 voxels.

Exercises

Levels: ★ apply the course, ★★ combine results or prove a short statement, ★★★ open-ended or research-flavoured. Exercises 6, 11, 14, 15, 17, 18 and 19 are computational.

1. ★ A snake shrinks exponentially. A closed snake with \alpha = 0.1, \beta = 0 and a constant potential starts as the circle of radius 20, parametrised at constant speed on [0, 1]. Compute its radius at t = 1 and the time at which the radius is halved. Compare with the curve shortening flow of the same circle.

   Solution: by Proposition 1.7, r(t) = 20 \, e^{-4\pi^2 \cdot 0.1 \, t}, so r(1) = 20 \, e^{-3.948} \approx 0.39. The radius is halved at t = \ln 2 / 3.948 \approx 0.176, whatever the initial radius.

   Under curve shortening, r^2 = 400 - 2t: the circle vanishes at t = 200 and halves at t = 150. The parametric snake shrinks at a rate independent of the size; the geometric flow shrinks small circles faster than large ones.

2. ★ Choosing the time step. For a snake of n = 100 nodes with a = 0.1 and b = 0.05, give the largest stable explicit time step, the smallest nonzero eigenvalue of A, and the time step allowed by the semi-implicit scheme.

   Solution: \lambda_{\max} = 4a + 16b = 1.2, so the explicit scheme needs \tau \leq 2/1.2 \approx 1.67 (Theorem 1.5).

   With \omega_1 = 4 \sin^2(\pi/100) = 0.003947, the smallest nonzero eigenvalue is \lambda_1 = a \omega_1 + b \omega_1^2 \approx 3.95 \times 10^{-4} . The semi-implicit scheme is stable for every \tau; its step is limited only by the external force, which should move a node by at most about one pixel.

3. ★ Finite elements with four nodes. For a closed curve of length 1 with n = 4 P1 elements (h = 1/4), write M and K (elasticity only, \alpha = 1) and compute the generalised eigenvalues of the lowest nonzero mode for finite elements and for finite differences.

   Solution: M = \frac{1}{24}\,\mathrm{tridiag}(1, 4, 1) and K = 4\,\mathrm{tridiag}(-1, 2, -1), cyclic. For k = 1, \theta = \pi/2 (Proposition 2.4): \lambda^{\mathrm{FD}} = 64 \sin^2 \frac{\pi}{4} = 32, \qquad \lambda^{\mathrm{FE}} = 96 \cdot \frac{1 - 0}{2 + 0} = 48 . The exact value 4\pi^2 \approx 39.5 lies between them, as Proposition 2.4 says.

4. ★ A circle under the geometric active contour. With g = 1 and a balloon \nu = 0.1, the geometric active contour is \phi_t = (\kappa + \nu) |\nabla \phi|. Start from the circle of radius 10. When does it disappear?

   Solution: with \phi = |x| - r(t) (Example 3.4), -\dot r = 1/r + \nu. Separating variables, T = \int_0^{10} \frac{r \, dr}{1 + \nu r} = \frac{1}{\nu} \Big( 10 - \frac{\ln(1 + 10\nu)}{\nu} \Big) = 10 \, (10 - 10 \ln 2) \approx 30.7 , instead of 50 without the balloon.

5. ★ Fast Marching updates by hand. With h = P = 1, compute the update of Proposition 6.7 for (a, b) = (0, 0), (2, 2.5) and (1, 3). Then extend the reasoning to 3D, where the equation is \sum_{i=1}^3 \max(u - a_i, 0)^2 = 1, with three axis neighbours a_i = 0.

   Solution: (0 + \sqrt 2)/2 \approx 0.707; (4.5 + \sqrt{2 - 0.25})/2 \approx 2.911; and since |1 - 3| \geq 1, the one-sided value 2.

   In 3D with all a_i = 0: one active term would give u = 1, larger than the other neighbours, which are then active too; two terms give 1/\sqrt 2, same objection. All three are active, 3u^2 = 1 and u = 1/\sqrt 3 \approx 0.577 (Proposition 7.3 states the general rule).

6. ★ A 32-neighbour stencil. Add the moves (\pm 1, \pm 3), (\pm 3, \pm 1), (\pm 2, \pm 3), (\pm 3, \pm 2) to the 16-neighbour stencil of Code 5.11. Predict the maximal relative error with Proposition 5.10, then measure it on a 201 \times 201 grid, for points at distance more than 40 from the source.

   Hints and expected result: in the first octant the directions have angles 0°, 18.4°, 26.6°, 33.7° and 45°. The largest gap is \arctan(1/3) = 18.4°, so the error is \frac{1}{\cos(9.2°)} - 1 \approx 0.0131 . The measurement gives 0.0131, against 0.0275 for 16 neighbours.

7. ★ Chan–Vese by hand. A 1D signal is I = (0.1, 0.2, 0.9, 0.6, 0.8, 0.3), with pixels 3 to 5 inside. Compute c_1, c_2 and the two data terms (\lambda_1 = \lambda_2 = 1). Should pixel 4 or pixel 6 change side?

   Solution: c_1 = 0.767 and c_2 = 0.2 (Proposition 4.3). The data terms are 0.0467 inside and 0.02 outside.

   A pixel prefers the inside when I is closer to c_1 than to c_2, that is above (c_1 + c_2)/2 = 0.483. Pixel 4 (0.6) stays inside and pixel 6 (0.3) stays outside: the partition is stable without a length term.

8. ★ Refraction of a minimal path. Let P = 1 for y > 0 and P = 2 for y < 0. Find the minimal path from (0, 1) to (3, -1) and its weighted length, and compare with the straight segment.

   Solution: the path is straight in each half-plane and crosses the axis at a point (x, 0) that minimises f(x) = \sqrt{x^2 + 1} + 2 \sqrt{(3 - x)^2 + 1} . Setting f'(x) = 0 gives Snell’s law \sin\theta_1 = 2 \sin\theta_2, with angles to the normal. Numerically x \approx 2.477 and f \approx 4.928, against 5.408 for the segment. The path saves time by staying longer in the fast medium, as light does: the eikonal equation is the equation of geometrical optics.

9. ★★ Stationary circles of a geodesic active contour. Let g(x) = 1 + (|x| - 10)^2, an edge indicator whose valley is the circle of radius 10. Find the circles centred at the origin that are stationary for the geodesic active contour without balloon, and say which are stable.

   Solution: for a circle of radius r, \kappa = 1/r and \langle \nabla g, N \rangle = g'(r), so the normal speed is -(g/r + g'), which vanishes when (r g(r))' = 0. This is the condition that L_g = 2\pi r g(r) is stationary, as Theorem 3.6 predicts.

   It reads 3r^2 - 40 r + 101 = 0, with roots r \approx 9.95, a minimum of L_g (stable), and r \approx 3.38, a maximum (unstable). The stable contour sits slightly inside the valley: the length term pulls it inwards.

10. ★★ Region competition separates textures with equal means. Model region 1 by \mathcal{N}(m_1, \sigma_1^2) and region 2 by \mathcal{N}(m_2, \sigma_2^2). Write the speed of Proposition 4.7. With m_1 = m_2 = 0, \sigma_1 = 1 and \sigma_2 = 2, and a flat boundary, which pixels push the boundary into region 2?

    Solution: up to the curvature term, v = \log \frac{\sigma_2}{\sigma_1} - \frac{(I - m_1)^2}{2 \sigma_1^2} + \frac{(I - m_2)^2}{2 \sigma_2^2} . Here v = \ln 2 - \frac{3}{8} I^2, positive when |I| < 1.36. Small values are better explained by the narrow model, large ones by the wide model. Chan–Vese, which compares only means, cannot separate these two regions.

11. ★★ A robust projection on the shape model. In Code 4.11, the occluded pixels pull the projection. Estimate the coefficients of the two modes by least squares on the pixels farther than 20 pixels from the centre of the bite only, including a free constant. Then vary the number of modes in the plain projection from 1 to 5.

    Hints and expected result: with the mask, the intersection over union with the true ellipse rises to 0.950, against 0.867 for the plain projection.

    The plain projection gives 0.735, 0.867, 0.855 and 0.828 with 1, 2, 3 and 5 modes. Beyond the two modes that the training set really has, the extra modes start fitting the occlusion: a shape prior must be both robust and low-dimensional.

12. ★★ Viscosity solutions at a ridge. Let p \neq q and U(x) = \min(|x - p|, |x - q|), the distance to two points. Show that U is a viscosity solution of |\nabla U| = 1 away from the two points, including on the bisector.

    Then show that V(x) = \max(|x - p|, |x - q|) satisfies the equation almost everywhere but is not a viscosity solution.

    Solution: off the bisector both functions are smooth with unit gradient. On the bisector, a C^1 function touching U from above has its gradient in the segment between the two unit gradients, whose norm is at most 1: U is a subsolution. No C^1 function touches the concave kink from below, so the supersolution test is empty.

    For V, the kink is convex. The function touching it from below with the average of the two unit gradients has a gradient of norm less than 1, so the supersolution test fails (compare Example 6.4).

13. ★★ Propagating from both ends. For a constant metric in 2D, compare the number of grid points made Alive when Fast Marching propagates from p until q is reached, at distance d = 100 pixels, and when it propagates from both points until the fronts meet. Redo the count in 3D, and explain why stopping when q becomes Alive is exact.

    Solution: about \pi d^2 \approx 31\,400 points against 2 \pi (d/2)^2 \approx 15\,700: half. In 3D, \frac43 \pi d^3 \approx 4.2 \times 10^6 against 1.05 \times 10^6: a quarter (Proposition 6.13).

    When q becomes Alive its value is final (Theorem 6.9), and back-propagation only needs the values below U(q), which are all Alive.

14. ★★ Fast Marching converges, Dijkstra does not. Compute the Euclidean distance to the centre of [-1, 1]^2 with the Fast Marching of Code 6.12 and with the 8-neighbour Dijkstra of Code 5.11, for n = 51, 101, 201, and report the maximal errors.

    Hints and expected result: Fast Marching gives 0.0437, 0.0263 and 0.0155; Dijkstra gives 0.0897, 0.0898 and 0.0898.

    The Dijkstra error does not decrease: it is the relative error 0.082 of Proposition 5.10 times the largest distance in the worst direction, about 1.08. Fast Marching is worse near the source and better everywhere else as soon as h is small.

15. ★★ Anisotropy and fixed stencils. Take the constant metric M = R \, \mathrm{diag}(1, \mu^2) \, R^\top, with R the rotation of 30°. Compute graph distances with 8 and 16 neighbours, each edge e costing \sqrt{e^\top M e}.

    Compare with the exact distance \sqrt{x^\top M x} for \mu = 1, 2, 5, 10, and justify the exact formula.

    Hints and expected result: for a constant metric, straight segments are minimal paths (the metric is a norm, and the triangle inequality of the norm shows that no broken line is shorter), so the exact distance is the norm of x.

    On an 81 \times 81 grid, for points farther than 20 from the source, the maximal relative errors are 0.082, 0.280, 1.135, 2.827 with 8 neighbours, and 0.027, 0.072, 0.221, 0.564 with 16 neighbours. Fixed stencils fail as the anisotropy grows (Note 7.8).

16. ★★ The factor 2 is attained. On the segment [0, 1] with P = 1, run farthest-point sampling from x_1 = 0 for k = 1 to 4 points, and compare the covering radii with the optimal ones.

    Solution: the samples are 0, 1, 0.5, then 0.25 (or 0.75), with covering radii 1, 0.5, 0.25 and 0.25. The optimal radii are 1/(2k): 0.5, 0.25, 0.167 and 0.125.

    The ratios are 2, 2, 1.5 and 2: the bound of Theorem 8.2 is attained, because the sampling starts at an end of the segment.

17. ★★★ Extracting a vascular tree. Build a synthetic Y-shaped vessel (a trunk that splits into two branches of different radii), with noise. Compute the minimal action from the root with a vesselness or flux metric, then back-propagate from the two leaves.

    Hints and expected result: the two paths must merge at the bifurcation and share the trunk exactly (Proposition 8.7); measure the length of the common part and the distance of the merging point to the true bifurcation.

    Then add the radius dimension of Code 7.11 and check that the radius drops at the bifurcation. Report how the result depends on \varepsilon and on the noise level. This is a natural start for a project on tubular trees.

18. ★★★ A path that goes straight through a crossing. Draw two lines crossing at 60°, with noise. Compare the isotropic minimal path between the two ends of one line with a path computed on a graph over (x, y, \theta) with 16 orientations, where a move goes along the current orientation and a change of orientation costs \beta.

    Hints and expected result: the isotropic path may switch lines at the crossing when the other line is cheaper. The lifted path keeps its orientation through the crossing when \beta is large enough. Find the smallest such \beta, and relate the construction to the curvature-penalised models of ch. 9.

19. ★★★ A global Chan–Vese minimiser. Implement the convex relaxation of Proposition 9.1 on the image of Code 4.8, with c_1 = 0.7 and c_2 = 0.3 fixed. Use a primal–dual algorithm such as that of Chambolle and Pock (2011).

    Hints and expected result: threshold u^\star at t = 0.3, 0.5, 0.7. The three regions should be almost identical, since u^\star is nearly binary, and their energy should not exceed that of the level-set result. Discuss where the level-set flow stopped in a local minimum, and what the relaxation cannot do (estimate c_1, c_2 jointly, include a shape prior).

Projects and assessment

The official assessment, per master-mva.com, is a project; the Université Paris-Saclay course detail classifies it as “Evaluation Continue non Intégrale” (partial continuous assessment), and no weights are displayed. The subjects below are sized for that project. Most of them study and reimplement a named research paper, in the MVA tradition, and every one has a measurable success criterion. The written exam and the rubric are proposals that show what a test of the course could look like; none of these proposals is official material.

Solo projects

1. [Solo · ≈ 20 h] Balloons and finite elements: a paper study. Study Cohen (1991) and Cohen and Cohen (1993), an official reference.
   • Implement: the balloon model with finite differences and with P1 finite elements, and two potentials, the gradient potential and a distance-map potential.
   • Data: synthetic images of disks, ellipses and a U-shaped object with a deep concavity, with Gaussian noise.
   • Success criterion: a mean distance below 1 pixel to the true boundary; a plot of the error against the number of nodes for both discretisations; capture-range maps for both potentials.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 1, 2.
2. [Solo · ≈ 20 h] Minimal paths with Fast Marching: a paper study. Reproduce the approach of Cohen and Kimmel (1997).
   • Implement: Fast Marching in 2D (Code 6.12 as a start), back-propagation by sub-pixel gradient descent with a second-order integrator, and partial propagation stopped at the target.
   • Data: synthetic curves with gaps and noise, and contours to close between two clicked points.
   • Success criterion: a centreline error below 1 pixel on synthetic curves; a comparison with Dijkstra on 4, 8 and 16 neighbours that shows the metrication error of Proposition 5.10.
   • Deliverables: code, report, talk. Chapters 5, 6.
3. [Solo · ≈ 20 h] Geodesic remeshing of an image domain: a paper study. Study Peyré and Cohen (2006).
   • Implement: farthest-point sampling with Fast Marching, the geodesic Voronoi diagram by front competition, and the dual Delaunay triangulation, for a metric built from the image gradient.
   • Data: a photograph from the scikit-learn sample images (load_sample_image("china.jpg")), converted to grey levels, or synthetic images with sharp edges.
   • Success criterion: the error of the piecewise-linear interpolation of the image on the mesh against the number of vertices, compared with a uniform mesh with the same number of vertices.
   • Deliverables: code, report, talk. Chapters 6, 8.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Vessels in space and radius: a paper study. Reproduce the idea of Li and Yezzi (2007) and the anisotropic metric of Benmansour and Cohen (2011), an official reference.
   • Data: synthetic 2D and 3D vessel trees with known centrelines and radii, then the public DRIVE retinal images (Staal and coauthors, 2004).
   • Work split: one member builds the multi-scale metrics (vesselness, oriented flux); one writes Fast Marching on the lifted grid, isotropic then anisotropic; one handles tree extraction from several leaves and the evaluation.
   • Success criterion: on synthetic data, centreline error below 1 pixel and radius error below 0.5 pixel; on DRIVE, extracted centrelines overlaid on the manual annotations, with the fraction of centreline points that fall on an annotated vessel pixel reported for each image (aim: above 95 %).
   • Deliverables: code, a 6-page report, a 15-minute talk. Chapters 6, 7, 8.
2. [Group of 2–3 · ≈ 30 h each] Region segmentation: local flows against global minimisers. Compare the level-set Chan–Vese flow (Chan and Vese, 2001) with the convex relaxation of Chan, Esedoglu and Nikolova (2006).
   • Data: synthetic images with weak edges, textures of equal means (Exercise 10) and occlusions, then grey-level images of the Berkeley segmentation data set (Martin and coauthors, 2001).
   • Work split: one member writes the level-set flows (geodesic active contour, Chan–Vese); one the primal–dual solver for the relaxation; one region competition with Gaussian models and a PCA shape prior.
   • Success criterion: a table of final energies and intersections over union, showing where the flows stop in a local minimum; the thresholded relaxation must never have a higher Chan–Vese energy than the flow, for the same fixed constants (Proposition 9.1).
   • Deliverables: code, report, talk. Chapters 3, 4, 9.
3. [Group of 2–3 · ≈ 30 h each] Virtual endoscopy path planning in 3D: a paper study. Study Deschamps and Cohen (2001).
   • Data: synthetic 3D tubular phantoms (a curved tube with folds and a varying radius), generated by the group, about 200^3 voxels.
   • Work split: one member writes 3D Fast Marching with partial and two-sided propagation; one designs the centring metric from the distance to the wall; one evaluates the paths and renders a fly-through as a sequence of slices.
   • Success criterion: a path whose distance to the wall stays above 80 % of the local radius, and a measured speed-up of two-sided propagation close to the factor 4 of Exercise 13.
   • Deliverables: code, report, talk. Chapters 6, 7, 8.

Written exam (proposal)

Format: 3 hours, one handwritten A4 sheet allowed, no electronic device; four exercises, independent of each other.

1. Course question (5 points). State and prove the level-set equation for a front moving with normal speed F. Then state the first variation of the weighted length \int g \, ds and deduce the level-set form of the geodesic active contour.

   Expected: Theorem 3.3 with the differentiation of \phi(C(p, t), t) = 0; the first variation of Theorem 3.6; then F = -(g\kappa + \langle \nabla g, N \rangle) substituted with N = \nabla\phi / |\nabla\phi|. Full marks need the role of the advection term explained.

2. Computation (5 points). Run Fast Marching by hand on a 3 \times 3 grid with h = 1, P = 1 and the source at the corner (0, 0). Give the order in which the nodes become Alive and the final values, and compare with the Euclidean distances and with a 4-neighbour Dijkstra.

   Expected: the two nodes adjacent to the source get 1, and the diagonal node (1, 1) gets 1.707. The two other corners on the axes get 2.

   Then the nodes (1, 2) and (2, 1) get 2.545, and the far corner gets 3.252.

   The exact values are 1.414, 2.236 and 2.828; Dijkstra gives 2, 3 and 4. Full marks need the one-sided and two-sided cases of Proposition 6.7 identified correctly.

3. Short proof (5 points). Let U be the action from two sources p and q. Prove that the minimum of U on the set where the two fronts meet is d_P(p, q)/2, attained at the midpoint of a minimal path. Explain how this speeds up the extraction of a minimal path between the two sources.

   Expected: the triangle inequality for the lower bound and the midpoint for the equality (Proposition 6.13); propagation from both ends stopped at the first meeting point, with half the area in 2D and a quarter of the volume in 3D for a constant metric.

4. Problem (5 points). A closed snake has n = 50 nodes, a = 0.2 and b = 0.1. (a) Give the spectrum of its matrix and the largest stable explicit step. (b) Without image force, how does a circle evolve under the continuous flow? (c) Write the balloon force and explain the choice k > k_1.

   Expected: (a) eigenvalues a\omega_k + b\omega_k^2 in [0, 2.4], so \tau \leq 0.83 (Theorem 1.5).

   2. The circle stays a circle and shrinks exponentially at the rate 4\pi^2\alpha + 16\pi^4\beta (Proposition 1.7). With h = 1/50, \alpha = a h^2 = 8 \times 10^{-5} and \beta = b h^4 = 1.6 \times 10^{-8}, the rate is 3.2 \times 10^{-3}, close to the smallest nonzero eigenvalue \lambda_1 of the matrix.

   3. Definition 1.8: a normal pressure balanced on edges by a normalised image force slightly stronger than the pressure.

Grading rubric (proposal)

References

Official references (master-mva.com)

• L. D. Cohen and I. Cohen. Finite element methods for active contour models and balloons for 2D and 3D images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 15(11), November 1993.
• L. D. Cohen. Minimal paths and fast marching methods for image analysis. In Mathematical Models in Computer Vision: The Handbook, Springer, 2005 (as listed on the course page; the handbook is edited by N. Paragios, Y. Chen and O. Faugeras).
• G. Peyré, M. Péchaud, R. Keriven and L. D. Cohen. Geodesic methods in computer vision and graphics. Foundations and Trends in Computer Graphics and Vision, 5(3–4):197–397, 2010.
• F. Benmansour and L. D. Cohen. Tubular structure segmentation based on minimal path method and anisotropic enhancement. International Journal of Computer Vision, 92(2):192–210, 2011.

Textbooks

• G. Aubert and P. Kornprobst. Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations. Springer, 2nd edition, 2006.
• M. Bardi and I. Capuzzo-Dolcetta. Optimal Control and Viscosity Solutions of Hamilton–Jacobi–Bellman Equations. Birkhäuser, 1997.
• P.-L. Lions. Generalized Solutions of Hamilton–Jacobi Equations. Pitman, 1982.
• S. Osher and R. Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer, 2003.
• J. A. Sethian. Level Set Methods and Fast Marching Methods. Cambridge University Press, 2nd edition, 1999.

Deformable models, level sets and regions (ch. 1–4)

• D. Adalsteinsson and J. A. Sethian. A fast level set method for propagating interfaces. Journal of Computational Physics, 118(2):269–277, 1995.
• V. Caselles, F. Catté, T. Coll and F. Dibos. A geometric model for active contours in image processing. Numerische Mathematik, 66(1):1–31, 1993.
• V. Caselles, R. Kimmel and G. Sapiro. Geodesic active contours. International Journal of Computer Vision, 22(1):61–79, 1997.
• T. F. Chan and L. A. Vese. Active contours without edges. IEEE Transactions on Image Processing, 10(2):266–277, 2001.
• L. D. Cohen. On active contour models and balloons. CVGIP: Image Understanding, 53(2):211–218, 1991.
• D. Cremers, M. Rousson and R. Deriche. A review of statistical approaches to level set segmentation: integrating color, texture, motion and shape. International Journal of Computer Vision, 72(2):195–215, 2007.
• M. Kass, A. Witkin and D. Terzopoulos. Snakes: active contour models. International Journal of Computer Vision, 1(4):321–331, 1988.
• M. E. Leventon, W. E. L. Grimson and O. Faugeras. Statistical shape influence in geodesic active contours. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2000.
• R. Malladi, J. A. Sethian and B. C. Vemuri. Shape modeling with front propagation: a level set approach. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(2):158–175, 1995.
• D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics, 42(5):577–685, 1989.
• S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.
• N. Paragios and R. Deriche. Geodesic active regions and level set methods for supervised texture segmentation. International Journal of Computer Vision, 46(3):223–247, 2002.
• A. Tsai, A. Yezzi, W. Wells, C. Tempany, D. Tucker, A. Fan, W. E. Grimson and A. Willsky. A shape-based approach to the segmentation of medical imagery using level sets. IEEE Transactions on Medical Imaging, 22(2):137–154, 2003.
• L. A. Vese and T. F. Chan. A multiphase level set framework for image segmentation using the Mumford and Shah model. International Journal of Computer Vision, 50(3):271–293, 2002.
• C. Xu and J. L. Prince. Snakes, shapes, and gradient vector flow. IEEE Transactions on Image Processing, 7(3):359–369, 1998.
• H.-K. Zhao, T. Chan, B. Merriman and S. Osher. A variational level set approach to multiphase motion. Journal of Computational Physics, 127(1):179–195, 1996.
• S. C. Zhu and A. Yuille. Region competition: unifying snakes, region growing, and Bayes/MDL for multiband image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 18(9):884–900, 1996.

Minimal paths, Fast Marching and metrics (ch. 5–7)

• G. Barles and P. E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equations. Asymptotic Analysis, 4(3):271–283, 1991.
• L. D. Cohen. Multiple contour finding and perceptual grouping using minimal paths. Journal of Mathematical Imaging and Vision, 14(3):225–236, 2001.
• L. D. Cohen and R. Kimmel. Global minimum for active contour models: a minimal path approach. International Journal of Computer Vision, 24(1):57–78, 1997.
• M. G. Crandall and P.-L. Lions. Viscosity solutions of Hamilton–Jacobi equations. Transactions of the American Mathematical Society, 277(1):1–42, 1983.
• T. Deschamps and L. D. Cohen. Fast extraction of minimal paths in 3D images and applications to virtual endoscopy. Medical Image Analysis, 5(4):281–299, 2001.
• A. F. Frangi, W. J. Niessen, K. L. Vincken and M. A. Viergever. Multiscale vessel enhancement filtering. Medical Image Computing and Computer-Assisted Intervention (MICCAI), Lecture Notes in Computer Science 1496, Springer, 1998.
• R. Kimmel and J. A. Sethian. Computing geodesic paths on manifolds. Proceedings of the National Academy of Sciences, 95(15):8431–8435, 1998.
• M. W. K. Law and A. C. S. Chung. Three dimensional curvilinear structure detection using optimally oriented flux. European Conference on Computer Vision (ECCV), Lecture Notes in Computer Science 5305, Springer, 2008.
• H. Li and A. Yezzi. Vessels as 4-D curves: global minimal 4-D paths to extract 3-D tubular surfaces and centerlines. IEEE Transactions on Medical Imaging, 26(9):1213–1223, 2007.
• J.-M. Mirebeau. Anisotropic fast-marching on Cartesian grids using lattice basis reduction. SIAM Journal on Numerical Analysis, 52(4):1573–1599, 2014.
• E. Rouy and A. Tourin. A viscosity solutions approach to shape-from-shading. SIAM Journal on Numerical Analysis, 29(3):867–884, 1992.
• J. A. Sethian. A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences, 93(4):1591–1595, 1996.
• J. A. Sethian and A. Vladimirsky. Ordered upwind methods for static Hamilton–Jacobi equations: theory and algorithms. SIAM Journal on Numerical Analysis, 41(1):325–363, 2003.
• J. N. Tsitsiklis. Efficient algorithms for globally optimal trajectories. IEEE Transactions on Automatic Control, 40(9):1528–1538, 1995.
• H. Zhao. A fast sweeping method for eikonal equations. Mathematics of Computation, 74(250):603–627, 2005.

Remeshing and applications (ch. 8)

• R. Ardon and L. D. Cohen. Fast constrained surface extraction by minimal paths. International Journal of Computer Vision, 69(1):127–136, 2006.
• F. Benmansour and L. D. Cohen. Fast object segmentation by growing minimal paths from a single point on 2D or 3D images. Journal of Mathematical Imaging and Vision, 33(2):209–221, 2009.
• T. F. Gonzalez. Clustering to minimize the maximum intercluster distance. Theoretical Computer Science, 38:293–306, 1985.
• V. Kaul, A. Yezzi and Y. Tsai. Detecting curves with unknown endpoints and arbitrary topology using minimal paths. IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(10):1952–1965, 2012.
• G. Peyré and L. D. Cohen. Geodesic remeshing using front propagation. International Journal of Computer Vision, 69(1):145–156, 2006.

Research directions (ch. 9)

• B. Appleton and H. Talbot. Globally minimal surfaces by continuous maximal flows. IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(1):106–118, 2006.
• E. J. Bekkers, R. Duits, A. Mashtakov and G. R. Sanguinetti. A PDE approach to data-driven sub-Riemannian geodesics in SE(2). SIAM Journal on Imaging Sciences, 8(4):2740–2770, 2015.
• Y. Boykov and V. Kolmogorov. Computing geodesics and minimal surfaces via graph cuts. IEEE International Conference on Computer Vision (ICCV), 2003.
• A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120–145, 2011.
• T. F. Chan, S. Esedoglu and M. Nikolova. Algorithms for finding global minimizers of image segmentation and denoising models. SIAM Journal on Applied Mathematics, 66(5):1632–1648, 2006.
• D. Chen, J.-M. Mirebeau and L. D. Cohen. Global minimum for a Finsler elastica minimal path approach. International Journal of Computer Vision, 122(3):458–483, 2017.
• K. Crane, C. Weischedel and M. Wardetzky. Geodesics in heat: a new approach to computing distance based on heat flow. ACM Transactions on Graphics, 32(5), 2013.
• H. Kervadec, J. Bouchtiba, C. Desrosiers, É. Granger, J. Dolz and I. Ben Ayed. Boundary loss for highly unbalanced segmentation. Medical Imaging with Deep Learning (MIDL), 2019.
• D. Marcos, D. Tuia, B. Kellenberger, L. Zhang, M. Bai, R. Liao and R. Urtasun. Learning deep structured active contours end-to-end. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2018.
• J.-M. Mirebeau. Fast-marching methods for curvature penalized shortest paths. Journal of Mathematical Imaging and Vision, 60(6):784–815, 2018.
• S. Peng, W. Jiang, H. Pi, X. Li, H. Bao and X. Zhou. Deep snake for real-time instance segmentation. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
• T. Pock, D. Cremers, H. Bischof and A. Chambolle. An algorithm for minimizing the Mumford–Shah functional. IEEE International Conference on Computer Vision (ICCV), 2009.
• S. Shit, J. C. Paetzold, A. Sekuboyina, I. Ezhov, A. Unger, A. Zhylka, J. P. W. Pluim, U. Bauer and B. H. Menze. clDice — a novel topology-preserving loss function for tubular structure segmentation. IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2021.
• S. R. S. Varadhan. On the behavior of the fundamental solution of the heat equation with variable coefficients. Communications on Pure and Applied Mathematics, 20(2):431–455, 1967.

Data sets for the projects

• D. Martin, C. Fowlkes, D. Tal and J. Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. IEEE International Conference on Computer Vision (ICCV), 2001.
• J. Staal, M. D. Abràmoff, M. Niemeijer, M. A. Viergever and B. van Ginneken. Ridge-based vessel segmentation in color images of the retina. IEEE Transactions on Medical Imaging, 23(4):501–509, 2004.