Level Set Methods

A level set method represents a moving contour implicitly, as the set of points where a function defined on the whole image grid equals zero, and moves the contour by evolving that function with a partial differential equation. Because the contour is never stored as a list of points, it can split into several pieces, merge with others, and develop corners without any bookkeeping. Osher and Sethian introduced the approach in 1988 for fronts in physics, such as flames and crystal growth; in image segmentation the geometric and geodesic active contours of the 1990s and region-based models such as Chan–Vese were all formulated and implemented with it. The Active Contours article covers parametric snakes, balloon forces and gradient vector flow; this article covers the implicit representation: the level set equation and how it is solved on a grid, geodesic active contours, the Chan–Vese model and its relation to Mumford–Shah, what topological freedom costs, and the morphological approximations in scikit-image.

Problem

A curve $C(t)$ in the plane moves along its normal with a speed $F$ that may depend on its curvature, on the image, or on statistics of the regions it separates. Track it on a fixed pixel grid, through any number of splits and merges, and stop when it reaches a boundary or a minimum of an energy.

Osher and Sethian (1988) describe the two older ways of doing this and why both fail for large motions. Marker points placed along the curve bunch up where curvature builds, so the curve must be regridded, and the regridding adds an error that resembles diffusion; when two regions merge, parts of the boundary must be removed by ad hoc rules. Volume-of-fluid methods, which store the fraction of each cell that lies inside, handle topology easily but make the curvature of the front hard to compute. The level set formulation keeps the fixed grid of the second and the geometric accuracy of the first.

Inputs and Outputs

The input is an image, an initial contour (one or more closed curves, a seed region, or a preset such as a checkerboard), and a speed law, either written down directly or derived as the gradient descent of an energy. Implementations store the initial contour as a function $\phi_0$ on the grid, usually the signed distance to the curve.

The output is the final function $\phi$. Its zero level set is the contour, at sub-pixel resolution if interpolated, and the set where $\phi > 0$ (in the convention used below) is a binary mask whose connected components are the segmented regions. The method does not label the components; counting or separating them is a connected-component pass afterwards.

Intuition

Picture the image as a flat sheet and a surface hovering above it, with height $\phi(x, y)$. The contour is the coastline where the surface meets the sheet. Instead of moving the coastline, raise and lower the surface: where the surface rises the land grows, where it sinks the sea advances. A single bump with two peaks whose saddle sinks below sea level becomes two islands; two islands whose surfaces rise between them become one. Nothing in the surface itself changes character at these moments, which is why topology changes come free.

The price is that a one-dimensional curve is now carried by a two-dimensional function (a surface in 3D by a volume), and the function must be kept well behaved away from the coastline, where nothing constrains it.

Mathematical Formulation

The level set equation

Let $\phi(\mathbf{x}, t)$ be positive inside the contour and negative outside, with $C(t) = \{\mathbf{x} : \phi(\mathbf{x}, t) = 0\}$. Sign conventions differ between sources: Chan and Vese and scikit-image use positive inside, Kichenassamy et al. and the signed distance function article negative inside, and Osher and Sethian’s original paper puts the front at $\phi = 1$ with $\phi > 1$ inside. A point $\mathbf{x}(t)$ on the front keeps $\phi(\mathbf{x}(t), t) = 0$, so by the chain rule $\phi_t + \nabla\phi \cdot \mathbf{x}' = 0$. Only the normal component of the velocity matters; with the outward normal $\mathbf{n} = -\nabla\phi / \lvert\nabla\phi\rvert$ and normal speed $F$,

$$\phi_t = F \, \lvert \nabla \phi \rvert ,$$

the form Adalsteinsson and Sethian write as $\phi_t - F\lvert\nabla\phi\rvert = 0$ with $\phi$ positive inside. Geometric quantities of every level curve come from $\phi$ directly. The curvature is

$$\kappa = \nabla \cdot \frac{\nabla \phi}{\lvert \nabla \phi \rvert} = \frac{\phi_{xx}\phi_y^2 - 2\phi_x\phi_y\phi_{xy} + \phi_{yy}\phi_x^2}{(\phi_x^2 + \phi_y^2)^{3/2}} ,$$

which with $\phi$ positive inside is negative on a convex contour, so the motion $\phi_t = \kappa\lvert\nabla\phi\rvert$ shrinks circles. Osher and Sethian derived the equation (their Eq. 2.13) by viewing the front as the set reached at time $t$, which turns the motion into a Hamilton–Jacobi equation with a second-order, parabolic right-hand side when $F$ depends on curvature; the same derivation applies in any number of dimensions with mean or Gaussian curvature.

Upwinding and the entropy condition

With constant speed, a front develops corners, and past a corner the equation has many weak solutions. Osher and Sethian adopt the entropy condition Sethian had proposed earlier for burning fronts, that a point once burnt stays burnt, which picks out the solution obtained as the limit of $F = 1 - \varepsilon K$ when $\varepsilon \to 0$, where Osher and Sethian’s curvature $K$ is positive on a convex front ($K = -\kappa$ in the convention here). Their tool is the link between Hamilton–Jacobi equations and hyperbolic conservation laws: monotone, upwind schemes from shock computation pick the entropy solution and do not oscillate near corners. Their simplest scheme (Eqs. 3.7 and 3.11), written here for $\phi$ positive inside and outward motion ($F > 0$), replaces $\lvert\nabla\phi\rvert$ by

$$\nabla^{+}\phi = \Big( \min(D^{-x}\phi, 0)^2 + \max(D^{+x}\phi, 0)^2 + \min(D^{-y}\phi, 0)^2 + \max(D^{+y}\phi, 0)^2 \Big)^{1/2},$$

where $D^{\pm x}$ are one-sided differences; for $F < 0$ the min and max swap. Each difference is taken from the side the front is coming from. Osher and Sethian record that the second author’s first attempt used the centred Lax–Friedrichs scheme, which created spurious waves at the boundary of the computational domain, and that this is what led them to upwind schemes. They also built second- and higher-order versions from essentially non-oscillatory (ENO) interpolation and measured roughly first- and second-order convergence on a test problem.

Curvature terms and the time step

Terms that depend on curvature are parabolic, not hyperbolic, and Osher and Sethian approximate all of them, including the $\lvert\nabla\phi\rvert$ that multiplies $\kappa$, by central differences (their Eq. 3.14): the curvature and the gradient must be computed the same way, because near the centre of concentric level curves both divide by a nearly vanishing $\lvert\nabla\phi\rvert$. A passive advection term $\mathbf{U} \cdot \nabla\phi$ is upwinded one coordinate at a time (Eq. 3.15). Time stepping is explicit (forward Euler or Runge–Kutta), so the step is limited by a CFL condition: for the scheme of Eq. 3.11 they give $1 \ge 2\,(\Delta t/\Delta x\,\lvert H_1\rvert + \Delta t/\Delta y\,\lvert H_2\rvert)$, where $H_1, H_2$ are the derivatives of the Hamiltonian $-F\lvert\nabla\phi\rvert$ with respect to $\phi_x, \phi_y$, so the front moves only a fraction of a cell per step, and they note that the parabolic curvature term makes the limit somewhat smaller. Their demonstrations include a seven-pointed star burning outward at unit speed, which develops corners that then open into a circle, and a wound spiral moving with $F = 1 - 0.01K$ that pinches off into an outer front, which expands, and an inner one, which collapses and vanishes, all without any code to detect the event.

Reinitialization

Only the zero level set has meaning, but the other level sets move too. Under image-driven speeds, neighbouring level curves are pushed at different rates, and $\phi$ becomes steep in some places and flat in others; both spoil the differences above and, for region models, the regularized delta function discussed below. The usual remedy is to replace $\phi$ from time to time by the signed distance to its current zero set, for which $\lvert\nabla\phi\rvert = 1$. Chan and Vese (CAM report, Eq. 13) do this, at every step, by solving the equation of Sussman, Smereka and Osher (1994) to steady state:

$$\psi_\tau = \operatorname{sign}(\phi)\,\big(1 - \lvert\nabla\psi\rvert\big), \qquad \psi(0, \cdot) = \phi .$$

Its zero set does not move, and away from it $\lvert\nabla\psi\rvert \to 1$; Chan and Vese found that a few iterations were enough near the front and give the Godunov-type upwind flux they used. Getreuer (2012) treats reinitialization as optional and notes that it does not change the segmentation boundary. The fast marching method, which the signed distance function article describes, is another way to rebuild the distance.

Narrow bands

A front of length $n$ cells in an $n \times n$ grid costs $O(n^2)$ per step, and $O(n^3)$ in three dimensions, although only cells near the front matter. Adalsteinsson and Sethian (1995; the details below are from their September 1994 preprint) restricted the update to a tube of points within a distance of the front, cutting the work to $O(nk)$ for a tube $k$ cells wide. The tube is built from the signed distance, computed only near the curve, and reused for as many steps as possible. Barrier points just inside the tube’s edges trigger a rebuild when one of them changes sign, meaning the front has come close to the edge. They compared three treatments of the tube’s edge and found that freezing the edge values was disastrous for curvature-dependent motion; the best was to compute derivatives only in the interior of the tube and extend the update to the edge points. A second motive was the extension problem: in some applications the speed is known only on the front, and a narrow band needs it only nearby.

Geodesic active contours

The geometric active contours of Caselles, Catté, Coll and Dibos (1993) and of Malladi, Sethian and Vemuri (1994, 1995) moved the zero level set with curvature plus a constant inflation, slowed by an edge-stopping function $g$, as Caselles, Kimmel and Sapiro (1997) describe them in their Eq. 15:

$$\phi_t = g(I)\,(c + \kappa)\,\lvert\nabla\phi\rvert, \qquad g = \frac{1}{1 + \lvert\nabla \hat{I}\rvert^p}, \; p = 1 \text{ or } 2,$$

with $\hat{I}$ a Gaussian-smoothed image. The curve stops only where $g$ is zero, that is, at an ideal edge. Written for $\phi$ positive inside, $c > 0$ inflates the contour; Caselles, Kimmel and Sapiro move each level curve by $c + \kappa$ along the inward normal, so their positive $c$ deflates (their function is negative inside, in our reading; flipping the sign of $\phi$ flips only the sign of $c$). Caselles, Kimmel and Sapiro showed that minimizing a snake energy with no rigidity term is equivalent to finding a curve of minimal weighted length $\int g \, ds$ (the derivation through Maupertuis’s principle is in the Active Contours article), and wrote the gradient descent of the weighted length on the level set (their Eqs. 14 and 18):

$$\phi_t = \lvert\nabla\phi\rvert \, \nabla\cdot\Big( g \frac{\nabla\phi}{\lvert\nabla\phi\rvert} \Big) + c\,g\,\lvert\nabla\phi\rvert = g\,(\kappa + c)\,\lvert\nabla\phi\rvert + \nabla g \cdot \nabla\phi .$$

The new term $\nabla g \cdot \nabla\phi$ is an advection along $\nabla g$, which points toward the middle of an edge from both sides, so it pulls the contour into the valley of $g$ instead of merely slowing it. The authors argue that this lets the contour stop on boundaries whose gradient varies strongly along their length, and across small gaps, and that it partly removes the need for the constant $c$, which they relate to Cohen’s balloon force. They report that their wrench example gave practically the same result with $c = 0$, only more slowly, but also that $c$ remains of importance when the curve starts inside an object; for the earlier model they note that $c$ is what lets a convex initial curve capture a non-convex object, since a convex curve stays convex under pure curvature motion. Kichenassamy, Kumar, Olver, Tannenbaum and Yezzi reached the same flow independently (Caselles et al., note 6): they replaced Euclidean arc length by $ds_\phi = \phi\, ds$ for an edge function $\phi$, computed the gradient flow that shrinks this modified length fastest, and obtained the level set equation with the extra term $\nabla\phi \cdot \nabla\Psi$ (their Eqs. 5–8, in the 1994 version read here). The 3D analogue moves a surface toward a minimal weighted area.

Caselles, Kimmel and Sapiro discretized their equation with central differences in space and forward differences in time, which they found stable for their equation, but they note that a large $c$ needs the more careful schemes of Osher and Sethian. Started as the image frame, their contour split around two wrenches a few pixels apart; started as two small circles inside two objects, the curves evolved independently and found the inner and outer boundaries of both.

Region-based models: Chan–Vese

Edge-based flows see the image only through $g$. Chan and Vese (CAM report of 1998, published 2001) proposed a model “without edges”, in which the stopping comes from how well two constants fit the image inside and outside the curve:

$$F(c_1, c_2, C) = \mu \,\mathrm{Length}(C) + \nu \,\mathrm{Area}(\mathrm{inside}(C)) + \lambda_1 \!\int_{\mathrm{inside}(C)}\! \lvert u_0 - c_1 \rvert^2 dx + \lambda_2 \!\int_{\mathrm{outside}(C)}\! \lvert u_0 - c_2 \rvert^2 dx .$$

In almost all their experiments $\lambda_1 = \lambda_2 = 1$ and $\nu = 0$; the report also allows the length to a power $p$. With the Heaviside function $H$ and its derivative $\delta$, the length is $\int \delta(\phi)\lvert\nabla\phi\rvert$, the area $\int H(\phi)$, and the fit terms are weighted by $H(\phi)$ and $1 - H(\phi)$. For fixed $\phi$ the optimal constants are the region means, $c_1 = \int u_0 H(\phi) / \int H(\phi)$ and $c_2$ likewise outside. For fixed constants, gradient descent in an artificial time gives (report, Eq. 12, $p = 1$)

$$\phi_t = \delta_\varepsilon(\phi) \Big[ \mu \, \nabla\cdot\frac{\nabla\phi}{\lvert\nabla\phi\rvert} - \nu - \lambda_1 (u_0 - c_1)^2 + \lambda_2 (u_0 - c_2)^2 \Big],$$

with a Neumann condition on the image border. The regularization matters. Chan and Vese compare a $\delta_{1,\varepsilon}$ supported on $[-\varepsilon, \varepsilon]$ with $\delta_{2,\varepsilon}(z) = \varepsilon / (\pi(\varepsilon^2 + z^2))$, the derivative of $H_{2,\varepsilon}(z) = \tfrac12(1 + \tfrac{2}{\pi}\arctan(z/\varepsilon))$, which is nonzero everywhere. With the compact one, the equation acts only on a few level curves near the front and they sometimes obtained only a local minimum; with the second, which they used for all reported results, every level curve feels the image, interior contours appear without a second initial curve, and the algorithm in their words tends to find a global minimizer. Fig. 3 of the report shows a hole opening inside a noisy torus from a single initial curve. Because the energy is nonconvex, they do not claim more.

Numerically, Chan and Vese discretize the divergence semi-implicitly, solve the linear system iteratively (one iteration per step in their results, five for noisy images), and reinitialize $\phi$ to a signed distance at each step. Getreuer (2012) follows the same semi-implicit scheme with a single Gauss–Seidel sweep per step, with the curvature coefficients regularized by $\eta = 10^{-8}$, a tolerance on the change in $\phi$ as the stopping rule, and the checkerboard $\phi_0 = \sin(\pi x/5)\sin(\pi y/5)$ as an initialization that converges fast, whereas a low-curvature initial curve, such as a large ellipse, can need thousands of iterations. With $\nu = 0$ and the checkerboard, which phase is called inside is arbitrary. Getreuer calls $\mu$ the most important parameter: in Chan and Vese’s examples, $\mu = 0$ detects every small object, and a large $\mu$ groups nearby objects into one region, which they relate to Kanizsa’s proximity grouping. Chan, Sandberg and Vese extended the model to vector-valued images by replacing $\lvert u_0 - c_i\rvert$ with a Euclidean norm over channels (Getreuer, §4).

Relation to Mumford–Shah

Mumford and Shah (1989) posed segmentation as finding a piecewise smooth approximation $f$ of the image $g$ with a discontinuity set $\Gamma$:

$$E(f, \Gamma) = \mu^2 \iint_R (f - g)^2 \, dx\, dy + \iint_{R - \Gamma} \lVert\nabla f\rVert^2 \, dx\, dy + \nu \lvert\Gamma\rvert ,$$

a “cartoon” of the image in their word, with smooth regions and sharp edges. Dropping any term makes the infimum zero. They conjectured but could not prove that $E$ has a minimizer. Restricting $f$ to a constant on each region gives

$$E_0(\Gamma) = \sum_i \iint_{R_i} \big(g - \mathrm{mean}_{R_i} g\big)^2 \, dx\, dy + \nu_0 \lvert\Gamma\rvert ,$$

the limit of $E$ as $\mu \to 0$, whose minimization they showed to be well posed with geometric measure theory. For the full functional $E$, assuming a minimizer whose singular points are joined by $C^2$ arcs, they show that the only possible singularities are triple points where three arcs meet at 120°, crack tips, and arcs meeting the image border at right angles. The first variation of $E_0$ gives the equation $(\bar g^{+} - \bar g^{-})(\bar g^{+} + \bar g^{-} - 2g) + \nu_0 \,\mathrm{curv}(\Gamma) = 0$ along $\Gamma$, with $\bar g^{\pm}$ the means on either side.

Chan and Vese state that their model with $p = 1$, $\nu = 0$ and $\lambda_1 = \lambda_2$ is a particular case of this minimal partition problem, in which $f$ may take only two values and the boundary is the zero level set of one function. Expanding $(u_0 - c_2)^2 - (u_0 - c_1)^2 = (c_1 - c_2)(2u_0 - c_1 - c_2)$ shows that, on the front, the Chan–Vese steady state is the Mumford–Shah first-variation equation with $\bar g^{\pm} = c_1, c_2$ (our algebra; the sign of the curvature term depends on the orientation). The difference is in what is allowed: $E_0$ gives every region its own mean, while Chan–Vese gives all components of the inside one mean and all components of the outside another, so two objects of different brightness cannot both be fitted well, and, since every boundary point separates $\phi > 0$ from $\phi < 0$, triple junctions cannot occur (our reasoning). Getreuer illustrates the gap between the full piecewise smooth model and the two-phase one, and notes that the existence proof for Chan–Vese follows from Morel and Solimini’s proof for the piecewise constant case; Getreuer cites Vese and Chan’s later multiphase level set framework (2002) for more than two phases.

Algorithm

A generic level set segmentation:

1. Choose the speed: an explicit law such as the geodesic flow, or the gradient descent of an energy such as Chan–Vese.
2. Initialize $\phi$ as the signed distance to one or more initial curves (positive inside), or as a preset (checkerboard, disk).
3. Repeat until a stopping rule holds (change in $\phi$ below a tolerance, an iteration limit, or front motion below a threshold in the band):
   1. compute the image-dependent quantities: $g$ and $\nabla g$ once, or the region means $c_1$, $c_2$ at every step;
   2. update $\phi$ with one explicit or semi-implicit time step: upwind differences for constant-speed and advection terms, central differences for curvature, a time step within the CFL limit; with a narrow band, only cells in the band are updated;
   3. every few steps, reinitialize $\phi$ to a signed distance and, with a narrow band, rebuild the band.
4. Return the mask $\phi > 0$ and, if needed, the zero level set as sub-pixel polygons or, in 3D, a triangle mesh extracted by marching cubes.

Parameters

• Speed constant $c$ (geodesic and geometric models). Pushes the contour through flat regions and lets a convex initial curve capture non-convex objects. Caselles, Kimmel and Sapiro write, of the earlier geometric model, that it must be fixed in advance and that choosing it is not trivial, and that their gradient term reduces the need for it. Too large, and it carries the contour across weak edges; the example below shows both failures.
• Edge function $g$. The smoothing scale $\sigma$ of $\hat I$, the exponent $p$, or in scikit-image’s inverse_gaussian_gradient the steepness alpha (default 100) and sigma (default 5) in $g = 1/\sqrt{1 + \alpha \lvert\nabla G_\sigma * I\rvert}$. In our reading, a larger $\sigma$ widens the valley around an edge but also lowers $g$ in narrow gaps between objects.
• Length weight $\mu$ (Chan–Vese). Sets the scale of detected objects and the smoothness of the boundary. Chan and Vese note that its units depend on the image range, so a fixed value means different things for images scaled to 0–255 or 0–1; scikit-image rescales the image to $[0, 1]$ first.
• Fit weights $\lambda_1, \lambda_2$ and area weight $\nu$. Unequal $\lambda$ lets one phase have a wider spread of values (scikit-image documentation); Chan and Vese needed $\lambda_1 < \lambda_2$ for objects of different intensities. A nonzero $\nu$ makes the curve prefer shrinking or growing and is meaningful only when inside and outside are prescribed (Getreuer, §5.3).
• Regularization $\varepsilon$, time step and stopping tolerance. $\varepsilon = 1$ is Getreuer’s and scikit-image’s choice. In a semi-implicit scheme the time step multiplies $\delta_\varepsilon$, so scaling one is equivalent to scaling the other (our reading of the update).
• Initialization. The checkerboard converges fast and finds many objects; a single curve around one object segments that object if the regularized delta and the reinitialization keep the evolution local.
• Reinitialization frequency and band width. Every step in Chan and Vese, optional in Getreuer; with a narrow band, the band must be wide enough that the front cannot leave it between rebuilds.

Complexity

Each step of a full-grid method costs $O(N)$ for $N$ grid points (Getreuer notes that the Chan–Vese iteration is linear in the number of pixels; Adalsteinsson and Sethian write $O(n^2)$ for an $n \times n$ grid and $O(n^3)$ for a volume). A narrow band of width $k$ reduces this to $O(nk)$ in 2D. Explicit schemes need many small steps, since the CFL condition ties the step to the grid spacing and the speed; the number of steps depends on how far the front travels and, for Chan–Vese, strongly on the initialization. Reinitialization with the Sussman equation adds a few sweeps over the grid or the band. Region models also recompute the two means each step, which costs one pass over the image. The morphological approximations below replace floating-point updates with binary erosions and dilations and need no time step.

Implementation

scikit-image 0.25.2 provides three functions (source read at that tag):

• chan_vese implements the semi-implicit update of Getreuer’s Eq. 22 in a vectorized form: all neighbours are taken from the previous iterate, a Jacobi rather than Gauss–Seidel sweep, which the docstring explains as a speed trade-off in Python, with all four curvature coefficients evaluated at the centre pixel (a notation the docstring credits to Rami Cohen). It omits the area term $\nu$, computes $c_1, c_2$ with the hard threshold $\phi > 0$, uses $\delta(x) = \varepsilon/(\varepsilon^2 + x^2)$ without Getreuer’s factor $1/\pi$, does not reinitialize (the reinitialization function is a placeholder that returns $\phi$), rescales the image to $[0, 1]$, and stops when the root-mean-square change in $\phi$ falls below tol (default $10^{-3}$) or after max_num_iter (500). Defaults are mu=0.25, lambda1=lambda2=1, dt=0.5 and the checkerboard initialization; the 'disk' preset is a cone $(r - d)/r$ whose radius $r$ is the distance from the centre pixel to the nearest image edge (63 on a $128 \times 128$ image).
• morphological_chan_vese and morphological_geodesic_active_contour work on a binary array (described under Variants).

The example below implements Chan–Vese and a geodesic active contour on the level set in NumPy and compares them with the three library functions on two synthetic $128 \times 128$ images with Gaussian noise of standard deviation 0.2 on a contrast of 1: image A holds a disk and an annulus 10 pixels apart (two components and one hole), image B a dumbbell, two disks joined by a bar 13 pixels thick (one component). It tests whether the NumPy Chan–Vese reproduces scikit-image; how the regularized delta, the initialization and reinitialization change the result; whether the geodesic contour splits around the two objects of A, starting outside, and merges two seeds inside B, for several values of $c$; how far $\lvert\nabla\phi\rvert$ drifts from 1 without reinitialization; and whether a narrow band changes the result. Each mask is scored by its intersection over union with the truth, its number of 8-connected components and its number of holes. The geodesic contour uses central differences for curvature, the Osher–Sethian upwind gradient for the $c\,g\lvert\nabla\phi\rvert$ term and first-order upwinding for $\nabla g \cdot \nabla\phi$, with $g$ from inverse_gaussian_gradient(alpha=100, sigma=2), a time step of 0.2, 1,500 iterations, and 10 Sussman iterations (step 0.5, with the sign smoothed as $\phi/\sqrt{\phi^2 + 1}$, a choice of this example) every 20 steps. The narrow band restricts only the update to $\lvert\phi\rvert < 6$; reinitialization still runs on the whole grid, a simplification of Adalsteinsson and Sethian’s tube. The script runs in about 10 seconds with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, single-threaded (OMP_NUM_THREADS=1), and its output was identical in two runs. NumPy also prints warnings about empty means for the runs in which no front survives.

import numpy as np
from scipy.ndimage import distance_transform_edt as edt
from skimage.measure import label
from skimage.segmentation import (chan_vese, morphological_chan_vese,
                                  morphological_geodesic_active_contour,
                                  inverse_gaussian_gradient, disk_level_set)

N = 128
rr, cc = np.mgrid[:N, :N].astype(float)


def disk(cy, cx, r):
    return (rr - cy) ** 2 + (cc - cx) ** 2 <= r ** 2


# Image A: a disk and an annulus (the annulus has a hole); Image B: a dumbbell
truth_A = disk(64, 38, 20) | (disk(64, 90, 22) & ~disk(64, 90, 10))
truth_B = disk(64, 36, 18) | disk(64, 92, 18) | ((abs(rr - 64) <= 6) & (cc >= 36) & (cc <= 92))
img_A = truth_A + 0.2 * np.random.default_rng(0).standard_normal((N, N))
img_B = truth_B + 0.2 * np.random.default_rng(1).standard_normal((N, N))


def scores(mask, truth):
    iou = (mask & truth).sum() / (mask | truth).sum()
    comps = label(mask, connectivity=2).max()
    bg = label(~mask, connectivity=1)
    border = set(np.unique(np.r_[bg[0], bg[-1], bg[:, 0], bg[:, -1]]))
    holes = len(set(np.unique(bg)) - border - {0})
    return iou, comps, holes


def shifts(phi):
    P = np.pad(phi, 1, mode="edge")
    return P, P[1:-1, 2:], P[1:-1, :-2], P[2:, 1:-1], P[:-2, 1:-1]  # right, left, down, up


def sdf(inside):
    """Signed distance on the pixel grid, positive inside."""
    return edt(inside) - edt(~inside)


# ---------- Chan-Vese: semi-implicit (Jacobi) step of Getreuer Eq. 22, as in skimage ----------
def chan_vese_np(img, phi, mu=0.25, dt=0.5, iters=500, tol=1e-3, pi_delta=True,
                 reinit_every=0):
    f = (img - img.min()) / (img.max() - img.min())
    eps = 1.0
    delta = (lambda x: eps / (np.pi * (eps ** 2 + x ** 2))) if pi_delta else \
            (lambda x: eps / (eps ** 2 + x ** 2))
    for it in range(1, iters + 1):
        P, R, L, D, U = shifts(phi)
        cx, cy = (R - L) / 2, (D - U) / 2
        C1 = 1 / np.sqrt(1e-16 + (R - phi) ** 2 + cy ** 2)
        C2 = 1 / np.sqrt(1e-16 + (phi - L) ** 2 + cy ** 2)
        C3 = 1 / np.sqrt(1e-16 + cx ** 2 + (D - phi) ** 2)
        C4 = 1 / np.sqrt(1e-16 + cx ** 2 + (phi - U) ** 2)
        inside = phi > 0
        c1 = f[inside].mean() if inside.any() else 0.0
        c2 = f[~inside].mean() if (~inside).any() else 0.0
        d = dt * delta(phi)
        new = (phi + d * (mu * (R * C1 + L * C2 + D * C3 + U * C4)
                          - (f - c1) ** 2 + (f - c2) ** 2)) / (1 + mu * d * (C1 + C2 + C3 + C4))
        change = np.sqrt(((new - phi) ** 2).mean())
        phi = new
        if reinit_every and it % reinit_every == 0:
            phi = sdf(phi > 0)
        if tol and change < tol:
            break
    return phi > 0, it


# ---------- Geodesic active contour on the level set (CKS 1997, Eq. 18), phi > 0 inside ----------
def reinit(phi, steps=10, dtau=0.5):
    """Sussman-Smereka-Osher reinitialization, Godunov flux as in Chan & Vese (CAM 98-53, Eq. 13)."""
    s = phi / np.sqrt(phi ** 2 + 1.0)
    pos = phi > 0
    psi = phi.copy()
    for _ in range(steps):
        P, R, L, D, U = shifts(psi)
        a, b, c, d = psi - L, R - psi, psi - U, D - psi
        Gp = np.sqrt(np.maximum(np.maximum(a, 0) ** 2, np.minimum(b, 0) ** 2)
                     + np.maximum(np.maximum(c, 0) ** 2, np.minimum(d, 0) ** 2)) - 1
        Gm = np.sqrt(np.maximum(np.minimum(a, 0) ** 2, np.maximum(b, 0) ** 2)
                     + np.maximum(np.minimum(c, 0) ** 2, np.maximum(d, 0) ** 2)) - 1
        psi = psi - dtau * s * np.where(pos, Gp, Gm)
    return psi


def gac(g, phi, c, iters=1500, dt=0.2, reinit_every=20, band=None):
    gy, gx = np.gradient(g)
    updated = []
    active = np.ones_like(phi, bool)
    for it in range(iters):
        if band is not None and it % reinit_every == 0:
            active = np.abs(phi) < band
        P, R, L, D, U = shifts(phi)
        Dmx, Dpx, Dmy, Dpy = phi - L, R - phi, phi - U, D - phi
        px, py = (R - L) / 2, (D - U) / 2
        pxx, pyy = R - 2 * phi + L, D - 2 * phi + U
        pxy = (P[2:, 2:] - P[2:, :-2] - P[:-2, 2:] + P[:-2, :-2]) / 4
        curv = (pxx * py ** 2 - 2 * px * py * pxy + pyy * px ** 2) / (px ** 2 + py ** 2 + 1e-10)
        a = c * g                                    # balloon speed, outward if > 0
        grad_out = np.sqrt(np.minimum(Dmx, 0) ** 2 + np.maximum(Dpx, 0) ** 2
                           + np.minimum(Dmy, 0) ** 2 + np.maximum(Dpy, 0) ** 2)
        grad_in = np.sqrt(np.maximum(Dmx, 0) ** 2 + np.minimum(Dpx, 0) ** 2
                          + np.maximum(Dmy, 0) ** 2 + np.minimum(Dpy, 0) ** 2)
        balloon = np.maximum(a, 0) * grad_out + np.minimum(a, 0) * grad_in
        adv = gx * np.where(gx < 0, Dmx, Dpx) + gy * np.where(gy < 0, Dmy, Dpy)
        step = dt * (g * curv + adv + balloon)
        phi = phi + np.where(active, step, 0.0)
        updated.append(active.mean())
        if reinit_every and (it + 1) % reinit_every == 0:
            phi = reinit(phi)
    return phi, float(np.mean(updated))


def grad_dev(phi, width=3.0):
    gy, gx = np.gradient(phi)
    near = np.abs(phi) < width
    return float(np.abs(np.hypot(gx, gy)[near] - 1).mean())


def circle_phi(centres, r):
    return np.max([r - np.hypot(rr - cy, cc - cx) for cy, cx in centres], axis=0)


def iou(a, b):
    return (a & b).sum() / max((a | b).sum(), 1)


# ================= Chan-Vese and the morphological variants =================
rows = []
init_disk = (63 - np.hypot(rr - 63, cc - 63)) / 63                       # skimage 'disk' preset
init_chk = np.sin(np.pi / 5 * rr) * np.sin(np.pi / 5 * cc)               # Getreuer Eq. 15
for init, name in [("checkerboard", "checkerboard"), ("disk", "disk")]:
    ref, _, en = chan_vese(img_A, init_level_set=init, extended_output=True)
    rows.append(("A", f"skimage chan_vese, {name}", ref, truth_A, f"{len(en)} its"))
    phi0 = init_chk if init == "checkerboard" else init_disk
    for pi_delta, dn in [(False, "no 1/pi"), (True, "1/pi")]:
        m, it = chan_vese_np(img_A, phi0, pi_delta=pi_delta)
        rows.append(("A", f"NumPy CV, {name}, delta {dn}", m, truth_A,
                     f"{it} its, IoU vs skimage {iou(m, ref):.4f}"))
m = chan_vese(img_A, init_level_set="disk", tol=0)
rows.append(("A", "skimage chan_vese, disk, tol=0", m, truth_A, "500 its"))

step = np.where(np.hypot(rr - 64, cc - 38) <= 25, 1.0, -1.0)   # Getreuer Eq. 16, left disk only
for every, name in [(0, "CV left-disk init, no reinit"), (50, "CV left-disk init, SDF every 50")]:
    m, _ = chan_vese_np(img_A, step, iters=1000, tol=0, reinit_every=every)
    rows.append(("A", name, m, truth_A, "1000 its"))

gA = inverse_gaussian_gradient(img_A, alpha=100, sigma=2)
gB = inverse_gaussian_gradient(img_B, alpha=100, sigma=2)
big = circle_phi([(64, 64)], 58)
seeds = circle_phi([(64, 36), (64, 92)], 5)
rows.append(("A", "morphological_chan_vese", morphological_chan_vese(
    img_A, 200, init_level_set=disk_level_set(img_A.shape), smoothing=1).astype(bool), truth_A, ""))
rows.append(("A", "morph. GAC, balloon -1", morphological_geodesic_active_contour(
    gA, 400, big > 0, smoothing=1, balloon=-1).astype(bool), truth_A, ""))
rows.append(("B", "morph. GAC, balloon +1", morphological_geodesic_active_contour(
    gB, 400, seeds > 0, smoothing=1, balloon=1).astype(bool), truth_B, ""))

print(f"{'img':3s} {'method':36s} {'IoU':>6s} {'comps':>5s} {'holes':>5s}  note")
for img, name, m, t, note in rows:
    s = scores(m, t)
    print(f"{img:3s} {name:36s} {s[0]:6.3f} {s[1]:5d} {s[2]:5d}  {note}")

# ================= geodesic active contour level set =================
one = np.ones((N, N))                     # known answers with g = 1 (dt = 0.2)
r_eq = lambda phi: np.sqrt((phi > 0).sum() / np.pi)
print(f"\nscheme check, circle r=30: speed -1 for t=10 -> r {r_eq(gac(one, circle_phi([(64, 64)], 30), -1.0, iters=50)[0]):.2f}"
      f", speed +1 -> r {r_eq(gac(one, circle_phi([(64, 64)], 30), 1.0, iters=50)[0]):.2f}"
      f", curvature flow t=200 -> r {r_eq(gac(one, circle_phi([(64, 64)], 30), 0.0, iters=1000)[0]):.2f}"
      f" (exact {np.sqrt(900 - 400):.2f})")
print(f"\n{'GAC run':30s} {'c':>5s} {'IoU':>6s} {'comps':>5s} {'holes':>5s}  note")
for c in (-1.0, -0.5, -0.25, 0.0):
    phi_r, _ = gac(gA, big, c)
    phi_n, _ = gac(gA, big, c, reinit_every=0)
    phi_b, fb = gac(gA, big, c, band=6.0)
    for name, phi, note in [
            ("A split, reinit every 20", phi_r, f"||grad|-1| {grad_dev(phi_r):.3f}"),
            ("A split, no reinit", phi_n, f"||grad|-1| {grad_dev(phi_n):.3f}"),
            ("A split, narrow band 6", phi_b,
             f"updates {fb:.1%} of pixels, IoU vs full {iou(phi_b > 0, phi_r > 0):.4f}")]:
        s = scores(phi > 0, truth_A)
        print(f"{name:30s} {c:5.2f} {s[0]:6.3f} {s[1]:5d} {s[2]:5d}  {note}")
for c in (1.0, 0.5, 0.25):
    phi_m, _ = gac(gB, seeds, c)
    s = scores(phi_m > 0, truth_B)
    print(f"{'B merge, two seeds':30s} {c:5.2f} {s[0]:6.3f} {s[1]:5d} {s[2]:5d}")
edge = truth_A & ~(np.roll(truth_A, 1, 0) & np.roll(truth_A, -1, 0) & np.roll(truth_A, 1, 1) & np.roll(truth_A, -1, 1))
print(f"g on true boundary of A: mean {gA[edge].mean():.3f}; off the objects: 5th pct "
      f"{np.percentile(gA[~truth_A], 5):.3f}, median {np.median(gA[~truth_A]):.3f}")
print(f"truth A: 2 components, 1 hole; truth B: 1 component, 0 holes, {truth_B.mean():.1%} of the image")

img method                                  IoU comps holes  note
A   skimage chan_vese, checkerboard       0.939     2     1  58 its
A   NumPy CV, checkerboard, delta no 1/pi  0.939     2     1  58 its, IoU vs skimage 1.0000
A   NumPy CV, checkerboard, delta 1/pi    0.931     3     1  75 its, IoU vs skimage 0.9917
A   skimage chan_vese, disk               0.197     1     0  1 its
A   NumPy CV, disk, delta no 1/pi         0.197     1     0  1 its, IoU vs skimage 1.0000
A   NumPy CV, disk, delta 1/pi            0.197     1     0  1 its, IoU vs skimage 1.0000
A   skimage chan_vese, disk, tol=0        0.202     1     0  500 its
A   CV left-disk init, no reinit          0.509     1     0  1000 its
A   CV left-disk init, SDF every 50       0.510     1     0  1000 its
A   morphological_chan_vese               0.879     2     0
A   morph. GAC, balloon -1                0.846    16     0
B   morph. GAC, balloon +1                0.982     1     0

scheme check, circle r=30: speed -1 for t=10 -> r 19.55, speed +1 -> r 39.61, curvature flow t=200 -> r 21.71 (exact 22.36)

GAC run                            c    IoU comps holes  note
A split, reinit every 20       -1.00  0.000     0     0  ||grad|-1| nan
A split, no reinit             -1.00  0.000     0     0  ||grad|-1| nan
A split, narrow band 6         -1.00  0.000     0     0  updates 4.9% of pixels, IoU vs full 0.0000
A split, reinit every 20       -0.50  0.000     0     0  ||grad|-1| nan
A split, no reinit             -0.50  0.000     0     0  ||grad|-1| nan
A split, narrow band 6         -0.50  0.000     0     0  updates 9.5% of pixels, IoU vs full 0.0000
A split, reinit every 20       -0.25  0.494     1     0  ||grad|-1| 0.010
A split, no reinit             -0.25  0.167     1     0  ||grad|-1| 0.484
A split, narrow band 6         -0.25  0.494     1     0  updates 18.6% of pixels, IoU vs full 1.0000
A split, reinit every 20        0.00  0.248     1     0  ||grad|-1| 0.008
A split, no reinit              0.00  0.270     1     0  ||grad|-1| 0.666
A split, narrow band 6          0.00  0.248     1     0  updates 26.2% of pixels, IoU vs full 1.0000
B merge, two seeds              1.00  0.140     1     0
B merge, two seeds              0.50  0.140     1     0
B merge, two seeds              0.25  0.517     2     0
g on true boundary of A: mean 0.221; off the objects: 5th pct 0.291, median 0.646
truth A: 2 components, 1 hole; truth B: 1 component, 0 holes, 14.0% of the image

The NumPy Chan–Vese with scikit-image’s delta matches the library exactly on the checkerboard start: the same 58 iterations and an IoU of 1.0000 between the masks, which confirms the reading of the source above. Both find the two objects and the hole in the annulus. The checkerboard already places small contours everywhere, so this result does not test Chan and Vese’s claim that interior contours appear from a single initial curve; the disk runs were meant to. Adding Getreuer’s $1/\pi$ changes the run: 75 iterations, one extra small component, and 99.2% agreement with the library. Because $\delta$ enters the update only through the product $dt\,\delta(\phi)$, dividing it by $\pi$ is the same as dividing the time step by $\pi$; the two versions have the same fixed points, but the tolerance stops them at different moments (our reading of the update).

The disk preset shows how a stopping rule can masquerade as convergence. All three disk runs stop after one iteration with the initial disk as the answer (IoU 0.197). The cone $(63 - d)/63$ has $\lvert\nabla\phi\rvert = 1/63$, so the curvature coefficients $C_1, \ldots, C_4$ are about 63 each and the semi-implicit denominator $1 + \mu\,dt\,\delta\,(C_1 + \cdots + C_4)$ damps the first update below the tolerance (our reading of the code). With the tolerance switched off, 500 iterations move the IoU only to 0.202, consistent with Getreuer’s remark that low-curvature initial curves evolve slowly. Neither reinitialization run reproduces the effect Getreuer describes, a separate component appearing around a distant object: starting from a $\pm 1$ step function around the left disk, both runs keep only the left disk (IoU 0.509, the disk’s share of the truth). In our reading, where the step function is exactly flat the regularized curvature coefficients are $10^{8}$ and the semi-implicit update returns the average of the neighbours, so nothing can nucleate there; Getreuer’s example used a different image, $\mu$ and number of iterations.

morphological_chan_vese finds both objects but misses the hole (0 holes, IoU 0.879): it changes only pixels on the current boundary, where the gradient of the binary level set is nonzero, so a contour whose initial disk covers the annulus, hole included, could open the hole only by eroding through the ring (our reading of the code). The morphological geodesic contour merges the two seeds in the dumbbell (one component, IoU 0.982) but, deflating from the large circle, breaks into 16 pieces (IoU 0.846).

The scheme check comes close to the radii the equations predict for unit speed ($30 \pm 10$, less a small curvature contribution), and its curvature motion shrinks the circle about 3% too fast; errors of that size do not explain the geodesic contour’s failures below. With $c = -1$ or $-0.5$, the deflating contour vanishes. With $c = -0.25$ it ends as a single component (IoU 0.494) and does not split. With $c = 0$ the single component has an IoU of 0.248; if it still contains both objects, that is about 9,900 pixels, a circle of radius about 56 against the initial 58 (our arithmetic), so curvature alone shrinks the large circle only slowly. The merge test fails the other way: at $c = 0.5$ and 1, the inflating contour floods the image (IoU 0.140, the dumbbell’s share of the image), and at $c = 0.25$ the two seeds stay two separate components (IoU 0.517) and never meet. The printed values of $g$ suggest why: $g$ averages 0.221 on the true boundary, but noise pushes it down to 0.291 at the 5th percentile of the background, so, in our reading, a balloon strong enough to cross the low-$g$ pockets in the background is also strong enough to cross the edges. No value of $c$ in the grid produced the planned split or merge; Caselles, Kimmel and Sapiro’s own examples ran on real images with their own $g$ and stopping criteria.

Reinitialization keeps $\lvert\nabla\phi\rvert$ within about 1% of 1 near the front, while without it the mean deviation is fifty to eighty times larger; at $c = -0.25$ it also changes the outcome, so in this example reinitialization is part of the method, not only a numerical courtesy. The narrow band reproduces the full-grid masks exactly while updating a fifth to a quarter of the pixels per step; at $c = -1$ and $-0.5$ both versions are empty.

Several cases were added after the first run and are disclosed here. The first run used only the disk preset for Chan–Vese, and every disk run stopped after one iteration; the checkerboard initialization and the disk run without a tolerance were added for the second run. The first run tried the geodesic contour only with $c = -1$ (split) and $c = +1$ (merge), and both failed; the grid of $c$ values was then written into the plan before the second run, because the sources leave $c$ to the user, with the note “Grid fixed now, not to be changed after run 2”, and every value in it is printed above. The scheme check was added after the second run, following a diagnostic of the geodesic code on circles. The printed values of $g$, the dumbbell’s share of the image, the fixed iteration counts in the note column and the guard that prints an agreement of 0 when both masks are empty were added to the printout, which was also split into two tables, without changing any computed result.

Properties and Behavior

Topology changes and their cost

Splitting and merging need no special code: Osher and Sethian’s spiral pinches off, Caselles, Kimmel and Sapiro’s frame-shaped initial contour splits around two wrenches, and their two seeds separately find inner and outer boundaries. Caselles et al. stress that the number of objects need not be known in advance, unlike for parametric snakes, which need explicit procedures to split. The freedom has costs.

• No control over topology. The method cannot be told to return exactly one object or a curve without holes. Getreuer shows a separate component appearing around a second fish when the contour was meant to outline one, and uses reinitialization, which he notes may be needed to avoid including separate objects, to keep a single closed curve; with a global delta, Chan–Vese will also carve holes and detect objects far from the initial curve whenever the fit improves. When the topology is known in advance, this freedom is, in our reading, a liability, and a parametric snake, whose topology is fixed, may suit better.
• A front finds a hole only if it can reach it. A front that moves inward from outside, as in the geodesic and morphological runs above, cannot reach the inside of an annulus without crossing the ring; seeds started inside objects can wrap around a hole, as in Caselles, Kimmel and Sapiro’s outward example that found inner and outer boundaries; and Chan–Vese with a global delta can open a hole anywhere, because every level curve responds to the region fit.
• One more dimension of work. The curve is carried by a function on the whole grid, $O(n^2)$ per step in 2D and $O(n^3)$ in 3D, which is what narrow bands were invented to reduce.
• Upkeep of $\phi$. The embedding function drifts away from a distance function and must be reinitialized; reinitialization adds its own iterations.
• Unlabelled output. The components of $\phi > 0$ carry no identity; following a particular object through a split, or across video frames, needs a separate labelling or one level set function per object (our reasoning).

Other properties

Local minima. Geodesic flows and Chan–Vese descend nonconvex energies, and the result depends on the initialization; Getreuer shows that a specific initial curve can select a particular object. Caselles et al. interpret the term with $c$ as an area constraint with $c$ as a Lagrange multiplier, which they liken to minimizing the weighted length plus $c$ times the enclosed area.

Intrinsic and sub-pixel. The weighted length and its flow do not depend on how the curve is parameterized, and the zero level set can be interpolated between grid points; Caselles et al. note that this allows accurate measurements, such as a tumour’s area in their MR example.

Independence of the embedding. In the viscosity framework the evolution of the zero level set does not depend on which function embeds it (Caselles et al., citing Chen, Giga and Goto, and Evans and Spruck), which is why reinitialization is allowed at all.

Noise. Chan and Vese note that their model needs no smoothing of the image, even a very noisy one, so boundary locations are preserved; edge-based models, they point out, need a smoothed image for $g$, which also blurs the edges, and since discrete gradients are bounded, $g$ never reaches zero on an edge, so the curve may pass through.

Limitations

• Parameters interact. $c$, $g$’s scale, $\mu$, the time step, the reinitialization schedule and the stopping rule all change the result; the example shows a tolerance that stops on the first step and a balloon with no setting in the tested grid that both crosses noise and stops at edges.
• Edge-based models need edges. The geodesic contour stops only where $g$ is small, and a smoothed noisy image lowers $g$ in the background as well as on edges. Chan and Vese built their model for objects whose boundaries are not defined by gradient, such as blurred objects or groups of small objects; for textures whose means are the same inside and outside, they replaced the image by a feature image, such as the curvature or orientation of its level lines.
• Two phases with constant means. Chan–Vese assumes each phase is roughly constant; shading, texture or several objects of different brightness need multiphase, piecewise smooth or feature-based extensions.
• Numerical machinery. Upwinding, CFL limits, reinitialization and narrow-band bookkeeping are easy to get subtly wrong, and simple implementations, such as scikit-image’s Jacobi Chan–Vese without reinitialization, differ from the published algorithms in ways that change iteration counts and occasionally results.
• Speed. Explicit evolutions need many iterations; Getreuer reports thousands for low-curvature starts.

Against other segmentation methods, level sets find local optima of continuous energies from an initial guess, while graph cuts find a global optimum of a discrete two-label energy of the kind they can represent; Getreuer lists graph-cut algorithms for the Chan–Vese and Mumford–Shah energies (Zehiry et al.; Bae and Tai) among alternative solvers. The watershed partitions an image without iterations or parameters of this kind, and learned segmenters replace hand-designed speeds with features learned from labelled masks. Level sets remain useful where a smooth closed boundary, sub-pixel accuracy, 3D surfaces or topology changes matter, and as a way to refine a mask from another method.

Variants

Morphological snakes

Álvarez, Baumela and Márquez-Neila (CVPR 2010) and Márquez-Neila, Baumela and Álvarez (PAMI 2014) replaced the PDE by morphological operators on a binary level set. In scikit-image’s implementation, the level set $u$ is 0 or 1; curvature smoothing alternates two compositions of operators built from erosions and dilations with four 3-pixel line structuring elements in 2D (nine planar ones in 3D): a supremum of erosions followed by an infimum of dilations, and the reverse. In morphological_geodesic_active_contour, the balloon is a dilation or erosion by a $3 \times 3$ square applied where $g$ exceeds a threshold divided by $\lvert c\rvert$ (by default the threshold is the 40th percentile of $g$), and the attraction term sets $u$ to 1 or 0 according to the sign of $\nabla g \cdot \nabla u$. In morphological_chan_vese, pixels on the boundary of $u$ are set inside or outside according to the sign of $\lambda_1(I - c_1)^2 - \lambda_2(I - c_0)^2$. The scikit-image documentation states, citing the 2014 paper, that these operators are infinitesimally equivalent to the PDEs, need no time step and are faster; that paper could not be read for this article. The first author’s README advises trying the Chan–Vese variant first, because it works on the raw image and is more robust to noise, and the geodesic variant only for images with visible contours, after preprocessing such as inverse_gaussian_gradient. Because they change only boundary pixels, both behave like a compact delta: they do not open holes away from the front, as the example shows.

Extensions of the region model

Getreuer lists multiphase segmentation with several level set functions (Vese and Chan, 2002), nested curves, texture features from Gabor filters, and Tsai, Yezzi and Willsky’s curve-evolution framework for the full Mumford–Shah functional, as well as alternative solvers: a topological-derivative method (He and Osher), multigrid (Badshah and Chen) and graph cuts. Chan, Sandberg and Vese’s vector-valued version segments colour images jointly over channels, and the model is invariant to unitary transforms of the channels (Getreuer, §4).

Other uses of the representation

Sethian’s fast marching method, described in the signed distance function article, treats fronts that advance monotonically. TurboPixels grow superpixels from regularly spaced seeds by a level-set geometric flow. The same implicit representation underlies signed distance functions in 3D reconstruction, where the surface is again a zero level set. Caselles, Kimmel and Sapiro demonstrated the geodesic contour on an MR image of a tumour, a task of the kind medical image segmentation covers.