Active Contours

An active contour is a curve placed in an image that moves until it settles on a feature, usually an object boundary. Its motion minimizes an energy with two kinds of terms: internal terms that keep the curve smooth and image terms whose minima lie on edges, lines or other features. Kass, Witkin and Terzopoulos introduced the model in “Snakes: Active Contour Models” (International Journal of Computer Vision, 1988) and named the curves snakes after the way they slither while minimizing their energy. The term now covers a family: parametric snakes, which move a list of points; snakes with added pressure (balloons) or diffused image forces (gradient vector flow); and geometric contours, which move a curve defined implicitly and can split and merge. Active contours belong to the classical image segmentation methods and were used widely in medical imaging, where an operator places a rough outline and the snake finds the boundary (McInerney & Terzopoulos, 1996).

Kass et al. framed the snake as one way to keep low-level vision open to correction. An edge detector commits to edge pixels that must then be linked; a snake searches for a whole contour, and a user or a higher-level process can push it from one local minimum of its energy into another. The journal page is dated 1988, with a 1987 copyright line, and Terzopoulos’s publication list and Xu and Prince (1998) give the year as 1987; this article uses 1988, the year on the journal page.

Problem

Given an image $I(x, y)$ and a rough initial curve near an object, find a nearby curve that follows the object’s boundary, is smooth where the boundary is smooth, and bridges gaps where the boundary is missing or weak. Kass et al. did not try to find contours automatically. Their snake relies on something else, a user interface, an automatic attentional mechanism or a higher-level interpretation, to place it near the intended contour, and the energy minimization does the rest; they called this a “power assist” for a person pointing at a feature.

Inputs and Outputs

The input is an image, an initial curve given as an ordered list of points (closed for an object outline, open with fixed or free ends for a line), and the energy’s weights. The image enters through a scalar potential, or directly through a force field, computed once before the contour moves.

The output is the final list of points. For a closed contour, filling the polygon gives a region mask, but the snake itself is a boundary representation with sub-pixel coordinates, not a labelling of pixels.

Intuition

Think of the snake as a chain of beads joined by springs and stiffened like a thin rod. Kass et al. describe the first-order term as making the snake act like a membrane and the second-order term as making it act like a thin plate. The springs shorten the chain and the stiffness straightens it. The image is a landscape whose valleys lie on edges; each bead slides downhill. The curve stops where the pull of the landscape balances the pull of its neighbours.

Two consequences follow from this picture, and most of the method’s later history responds to them. First, a bead feels the landscape only where the landscape slopes: on a flat region far from any edge, the springs alone act and the chain shrinks, so the snake must start close to the boundary. Second, at the mouth of a narrow concavity the edges on either side pull the chain sideways, not inward, so the snake bridges the opening instead of entering it.

Mathematical Formulation

The energy

Kass et al. write the snake as a parametric curve $\mathbf{v}(s) = (x(s), y(s))$, $s \in [0, 1]$, with energy (their Eq. 1)

$$E^{*}_{\text{snake}} = \int_0^1 E_{\text{int}}(\mathbf{v}(s)) + E_{\text{image}}(\mathbf{v}(s)) + E_{\text{con}}(\mathbf{v}(s)) \, ds .$$

The internal energy (Eq. 2) is

$$E_{\text{int}} = \tfrac{1}{2}\left( \alpha(s) \lvert \mathbf{v}_s(s) \rvert^2 + \beta(s) \lvert \mathbf{v}_{ss}(s) \rvert^2 \right),$$

where subscripts denote derivatives with respect to $s$. The weight $\alpha$ penalizes stretching (tension) and $\beta$ penalizes bending (rigidity). Setting $\beta(s) = 0$ at a point lets the snake form a corner there. Kass et al. relate this “controlled continuity” spline to Tikhonov regularization.

The image energy (Eq. 3) combines three functionals:

$$E_{\text{image}} = w_{\text{line}} E_{\text{line}} + w_{\text{edge}} E_{\text{edge}} + w_{\text{term}} E_{\text{term}} .$$

• Lines. $E_{\text{line}} = I(x, y)$ (Eq. 4). The sign of $w_{\text{line}}$ decides whether the snake seeks light or dark lines.
• Edges. $E_{\text{edge}} = -\lvert \nabla I(x, y) \rvert^2$, so the snake is drawn to large image gradients. In practice the image is first smoothed with a Gaussian, $E = -\lvert \nabla (G_\sigma * I) \rvert^2$ (Xu & Prince, 1998, Eq. 3). Kass et al. also tried $-(G_\sigma * \nabla^2 I)^2$, whose minima lie on the zero crossings of the Marr–Hildreth detector (Eq. 5); the equation is labelled $E_{\text{line}}$ in the paper although the text calls it an edge functional.
• Terminations. $E_{\text{term}}$ is the curvature of the level lines of a slightly smoothed image $C = G_\sigma * I$ (Eqs. 6–8),

$$E_{\text{term}} = \frac{C_{yy} C_x^2 - 2 C_{xy} C_x C_y + C_{xx} C_y^2}{(C_x^2 + C_y^2)^{3/2}},$$

which attracts the snake to line ends and corners. Combined with the edge term, it let a snake settle on a subjective contour, the illusory outline in a standard Kanizsa-type figure, where the shape between the visible fragments is set by the smoothness term alone (their Fig. 5).

The constraint energy $E_{\text{con}}$ carries user forces. In the authors’ “Snake Pit” interface on a Symbolics Lisp Machine, a user could attach springs between a snake point and a fixed position, another snake or the mouse, and push snakes with a clipped $1/r^2$ repulsion drawn as a volcano icon.

The Euler–Lagrange equations

A local minimum of the energy satisfies, for constant $\alpha$ and $\beta$, the force balance written by Xu and Prince (1998, Eq. 6) as

$$\alpha \mathbf{v}_{ss} - \beta \mathbf{v}_{ssss} - \nabla E_{\text{ext}} = 0 ,$$

with $E_{\text{ext}} = E_{\text{image}} + E_{\text{con}}$. The first two terms are the internal force, the last the external force. Kass et al.’s appendix prints the continuous equations (10)–(11) with a plus sign on the $\alpha$ term; their discrete equation (14), $\alpha_i(\mathbf{v}_i - \mathbf{v}_{i-1}) - \alpha_{i+1}(\mathbf{v}_{i+1} - \mathbf{v}_i) + \ldots$, corresponds to the sign above, which is what the variational derivative of the energy gives.

Discretization and the semi-implicit step

With $n$ points $\mathbf{v}_i$ spaced $h$ apart, finite differences turn the internal energy into a quadratic form, and the Euler equations become (Kass et al., Eqs. 15–16)

$$A \mathbf{x} + \mathbf{f}_x(\mathbf{x}, \mathbf{y}) = 0, \qquad A \mathbf{y} + \mathbf{f}_y(\mathbf{x}, \mathbf{y}) = 0,$$

where $\mathbf{x}$ and $\mathbf{y}$ stack the coordinates, $\mathbf{f}_x = \partial E_{\text{ext}} / \partial x_i$, and $A$ is a pentadiagonal matrix (cyclic for a closed snake) built from $\alpha$ and $\beta$. To solve them, Kass et al. make the snake evolve in time and treat the two kinds of force differently: the internal forces, which are linear in the positions, implicitly, and the external forces, which depend nonlinearly on position, explicitly (Eqs. 17–20):

$$\mathbf{x}_t = (A + \gamma I)^{-1} \left( \gamma \mathbf{x}_{t-1} - \mathbf{f}_x(\mathbf{x}_{t-1}, \mathbf{y}_{t-1}) \right),$$

and the same for $\mathbf{y}$. At equilibrium the positions stop changing and the original equations hold. The paper calls $\gamma$ a step size; because it multiplies $\mathbf{x}_t - \mathbf{x}_{t-1}$, a larger $\gamma$ gives smaller moves, so it acts as an inverse time step (our reading of Eq. 17). Cohen (1991, Eq. 3) writes the same scheme with a time step $\tau$, which plays the role of $1/\gamma$, and the force $\mathbf{F} = -\nabla P$:

$$(I + \tau A)\, \mathbf{v}_t = \mathbf{v}_{t-1} + \tau \mathbf{F}(\mathbf{v}_{t-1}).$$

Why implicit? Kass et al. note that with a fully explicit scheme an impulse needs $O(n^2)$ iterations, each $O(n)$, to travel the length of the snake, so explicit snakes are “flaccid”; the semi-implicit scheme moves forces along the whole snake in one $O(n)$ iteration and tolerates rigid snakes and large steps. The price is that strong external forces still need small steps, because they are integrated explicitly.

Algorithm

1. Compute the external potential, for edges $P = -\lvert \nabla (G_\sigma * I) \rvert^2$, and its force field $-\nabla P$ on the pixel grid.
2. Initialize the snake as $n$ points near the target contour.
3. Build $A$ from $\alpha$, $\beta$ and $h$; factor $A + \gamma I$ once (an LU decomposition of a banded matrix takes $O(n)$).
4. Repeat until the points stop moving or an iteration limit is reached:
   1. interpolate the force at each point’s sub-pixel position;
   2. solve for the new positions with the factored matrix;
   3. optionally resample the points to keep their spacing even, and refactor.
5. Return the points.

Interpolating the force matters. Cohen (1991) observed that with forces known only at pixel centres a point can oscillate between the two pixels around a minimum, and used bilinear interpolation so that equilibria fall at zeros of a continuous force. Resampling is not in the original paper; Xu and Prince (1998) reparameterized their snakes dynamically to keep neighbouring points 0.5–1.5 pixels apart.

Parameters

• $\alpha$ (tension). Larger values make the snake shorter and smoother; with no image force the snake shrinks. scikit-image describes its alpha as making the snake contract faster.
• $\beta$ (rigidity). Larger values resist bending and round off corners; $\beta = 0$ allows corners. Xu and Prince used $\alpha = 0.6$, $\beta = 0$ in all their experiments and attribute the rounding of corners in their results to $\alpha$ alone.
• Image weights. $w_{\text{line}}$, $w_{\text{edge}}$ and $w_{\text{term}}$ select the features. Their scale relative to $\alpha$ and $\beta$ decides whether the image or the smoothness wins. Cohen reported good results with the internal weights of the order of $h^2$ for $\alpha$ and $h^4$ for $\beta$, so that all coefficients of the matrix have similar magnitude.
• $\sigma$ of the Gaussian. Larger $\sigma$ spreads the potential farther and enlarges the capture range but blurs and displaces boundaries and can erase a concavity (Xu & Prince, 1998).
• Time step ($\gamma$ or $\tau$). Too large and points jump across the edge they should stop on; too small and weak forces have no visible effect. Cohen noted that Kass et al. avoided the first problem by tuning the step by hand and replaced the tuning with a force of fixed magnitude (see Variants).
• Number of points and stopping rule. The number of points limits the detail the snake can represent. scikit-image’s documentation warns that its output has as many points as its input, so the initial snake needs enough points for the final contour.

Complexity

One iteration costs $O(n)$ for $n$ points with a banded factorization (Kass et al.), plus $O(n)$ force interpolations. Precomputing the potential costs one smoothing and one gradient over the image, $O(N)$ for $N$ pixels with a separable Gaussian. The number of iterations depends on how far the snake must travel and on the step size; there is no bound in the original paper, and practical implementations stop when the points stop moving. Gradient vector flow adds an iterative diffusion over the whole image, $O(N)$ per iteration; Xu and Prince ran about $\sqrt{N}$ iterations and reported, on 1998 hardware, MATLAB times of minutes for a $256 \times 256$ image.

Implementation

scikit-image implements the Kass snake as skimage.segmentation.active_contour (version 0.25.2, source read at that tag). It differs from the paper in several details worth knowing before using it:

• The image term is w_line * I + w_edge * sobel(I), the Sobel gradient magnitude, not its square, and the snake moves up the gradient of this term. Positive w_line therefore attracts the snake to bright regions, as the documentation says, which is the opposite sign convention to Kass et al.’s $E_{\text{line}} = I$ (our reading of the code).
• The term is interpolated with a quadratic bivariate spline, and $(A + \gamma I)$ is inverted once.
• Each step is capped at max_px_move pixels through a $\tanh$, which the paper does not have.
• The loop stops when the snake’s largest displacement, compared with each of its last ten positions, falls below convergence (default 0.1) in at least one comparison, or after max_num_iter (default 2500) iterations. Ends can be periodic, fixed or free.

The example below implements a snake from scratch in NumPy with Cohen’s form of the semi-implicit step and compares it with scikit-image on a synthetic image: a disk of radius 36 with a notch 16 pixels wide cut from the top, a boundary concavity of the kind Xu and Prince used, with Gaussian noise of standard deviation 0.05. It tests four claims from the sources: a traditional snake bridges a concavity; started inside a flat region with no pressure it shrinks; a balloon force lets it start far from the boundary; and gradient vector flow both extends the capture range and pulls the snake into the concavity. Each run is scored by its intersection over union with the true mask and by the fraction of the notch’s pixels that the snake leaves out (1 means the snake followed the notch; 0 means it bridged it). The snake uses Xu and Prince’s $\alpha = 0.6$, $\beta = 0$, a unit time step, edge and gradient vector flow forces scaled to a maximum magnitude of one, resampling to 1-pixel spacing every 5 iterations, and a fixed 1,500 iterations. Gradient vector flow follows their Eq. 16 with $\mu = 0.2$. It runs in about 6 seconds with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, single-threaded (OMP_NUM_THREADS=1).

import numpy as np
from scipy import ndimage as ndi
from scipy import sparse
from scipy.sparse.linalg import splu
from skimage.draw import polygon2mask
from skimage.filters import gaussian
from skimage.segmentation import active_contour

# Synthetic object: a disk with a deep notch (a boundary concavity)
H = W = 128
rr, cc = np.mgrid[:H, :W]
disk = (rr - 68) ** 2 + (cc - 64) ** 2 <= 36 ** 2
notch_box = (rr >= 25) & (rr <= 70) & (cc >= 56) & (cc <= 71)
mask = disk & ~notch_box
notch = disk & notch_box                      # pixels the snake should leave out
rng = np.random.default_rng(0)
img = mask.astype(float) + 0.05 * rng.standard_normal((H, W))
smooth = gaussian(img, sigma=2)

# Edge potential P = -|grad(G_sigma * I)|^2 and its negative gradient
gy, gx = np.gradient(smooth)
P = -(gx ** 2 + gy ** 2)
Py, Px = np.gradient(P)
mag = np.hypot(Px, Py)
trad = (-Px / mag.max(), -Py / mag.max())             # Xu-Prince: max magnitude 1
unit = (-Px / (mag + 1e-12), -Py / (mag + 1e-12))     # Cohen: normalised direction


def gvf(f, mu=0.2, iters=128, dt=1.0):
    """Gradient vector flow (Xu & Prince 1998, Eq. 16), unit pixel spacing."""
    fy, fx = np.gradient(f)
    b = fx ** 2 + fy ** 2
    u, v = fx.copy(), fy.copy()
    r = mu * dt
    for _ in range(iters):
        lu = ndi.laplace(u, mode="nearest")
        lv = ndi.laplace(v, mode="nearest")
        u = (1 - b * dt) * u + r * lu + b * fx * dt
        v = (1 - b * dt) * v + r * lv + b * fy * dt
    m = np.hypot(u, v).max()
    return u / m, v / m


edge_map = -P / (-P).max()                    # f = |grad(G*I)|^2 in [0, 1]
gvf128 = gvf(edge_map, iters=128)
gvf1000 = gvf(edge_map, iters=1000)


def circle(cy, cx, r, n):
    t = np.linspace(0, 2 * np.pi, n, endpoint=False)
    return np.stack([cy + r * np.sin(t), cx + r * np.cos(t)], axis=1)  # (row, col)


def resample(pts, spacing=1.0):
    closed = np.vstack([pts, pts[:1]])
    seg = np.hypot(*np.diff(closed, axis=0).T)
    s = np.concatenate([[0], np.cumsum(seg)])
    n = max(int(round(s[-1] / spacing)), 8)
    t = np.linspace(0, s[-1], n, endpoint=False)
    return np.stack([np.interp(t, s, closed[:, 0]), np.interp(t, s, closed[:, 1])], axis=1)


def snake_matrix(n, alpha, beta, tau):
    """Sparse LU of I + tau*A for a closed snake (A pentadiagonal, cyclic)."""
    e = sparse.identity(n, format="csc")
    d2 = sparse.diags([1, 1, -2, 1, 1], [-(n - 1), -1, 0, 1, n - 1], shape=(n, n), format="csc")
    A = -alpha * d2 + beta * (d2 @ d2)
    return splu((e + tau * A).tocsc())


def run_snake(init, force=None, balloon=0.0, k=0.5, alpha=0.6, beta=0.0, tau=1.0, iters=1500, edge=unit):
    """Semi-implicit snake: (I + tau A) v_t = v_{t-1} + tau F(v_{t-1})."""
    v = resample(init)
    for it in range(iters):
        if it % 5 == 0:
            v = resample(v)
            inv = snake_matrix(len(v), alpha, beta, tau)
        coords = [v[:, 0], v[:, 1]]
        if balloon:
            fx = ndi.map_coordinates(edge[0], coords, order=1, mode="nearest") * k
            fy = ndi.map_coordinates(edge[1], coords, order=1, mode="nearest") * k
            t = np.roll(v, -1, 0) - np.roll(v, 1, 0)                      # tangent (dy, dx)
            area = 0.5 * np.sum(v[:, 1] * np.roll(v[:, 0], -1) - np.roll(v[:, 1], -1) * v[:, 0])
            nrm = np.sign(area) * np.stack([-t[:, 1], t[:, 0]], 1)        # outward normal (dy, dx)
            nrm /= np.linalg.norm(nrm, axis=1, keepdims=True) + 1e-12
            F = np.stack([fy, fx], 1) + balloon * nrm
        else:
            fx = ndi.map_coordinates(force[0], coords, order=1, mode="nearest")
            fy = ndi.map_coordinates(force[1], coords, order=1, mode="nearest")
            F = np.stack([fy, fx], 1)
        v = inv.solve(v + tau * F)
        v = np.clip(v, 0, H - 1)
    return v


def score(snake):
    region = polygon2mask((H, W), snake)
    iou = (region & mask).sum() / (region | mask).sum()
    notch_rec = (notch & ~region).sum() / notch.sum()
    return iou, notch_rec


outside = circle(68, 64, 50, 200)
inside = circle(88, 64, 8, 50)

runs = {}
runs["skimage, outside"] = active_contour(smooth, outside, w_line=0, w_edge=1)
runs["skimage, no early stop"] = active_contour(smooth, outside, w_line=0, w_edge=1, convergence=0)
runs["edge force, outside"] = run_snake(outside, trad)
runs["edge force, inside"] = run_snake(inside, trad)
runs["inflate k1=0.25, inside"] = run_snake(inside, balloon=0.25)
runs["inflate k1=0.45, inside"] = run_snake(inside, balloon=0.45)
runs["deflate k1=0.25, outside"] = run_snake(outside, balloon=-0.25)
runs["deflate k1=0.45, outside"] = run_snake(outside, balloon=-0.45)
runs["inflate 0.25, max-norm, in"] = run_snake(inside, balloon=0.25, k=1.0, edge=trad)
runs["deflate 0.25, max-norm, out"] = run_snake(outside, balloon=-0.25, k=1.0, edge=trad)
runs["GVF 128 it, outside"] = run_snake(outside, gvf128)
runs["GVF 1000 it, outside"] = run_snake(outside, gvf1000)
runs["GVF 1000 it, inside"] = run_snake(inside, gvf1000)

print(f"{'run':28s} {'IoU':>6s} {'notch left out':>15s}")
for name, s in runs.items():
    iou, nr = score(s)
    print(f"{name:28s} {iou:6.3f} {nr:15.3f}")

print("\nforce magnitude below the object (bottom edge at row 104), column 64")
print(f"{'distance (px)':14s}" + "".join(f"{d:>7d}" for d in (2, 5, 10, 20)))
for name, (fx, fy) in [("edge force", trad), ("GVF 128 it", gvf128), ("GVF 1000 it", gvf1000)]:
    m = np.hypot(fx, fy)
    print(f"{name:14s}" + "".join(f"{m[104 + d, 64]:7.3f}" for d in (2, 5, 10, 20)))

run                             IoU  notch left out
skimage, outside              0.452           0.000
skimage, no early stop        0.833           0.064
edge force, outside           0.866           0.159
edge force, inside            0.000           1.000
inflate k1=0.25, inside       0.354           0.975
inflate k1=0.45, inside       0.183           0.432
deflate k1=0.25, outside      0.808           0.990
deflate k1=0.45, outside      0.455           0.974
inflate 0.25, max-norm, in    0.985           0.970
deflate 0.25, max-norm, out   0.971           0.990
GVF 128 it, outside           0.888           0.369
GVF 1000 it, outside          0.982           0.972
GVF 1000 it, inside           0.982           0.974

force magnitude below the object (bottom edge at row 104), column 64
distance (px)       2      5     10     20
edge force      0.662  0.012  0.000  0.001
GVF 128 it      0.746  0.498  0.195  0.010
GVF 1000 it     0.777  0.678  0.538  0.369

The force table is the capture range in numbers. Two pixels below the object, the edge force and both gradient vector flow fields are strong. Five pixels away the edge force has fallen to about 1% of its maximum, and at ten it is zero to three decimals: beyond the reach of the Gaussian, a snake point feels only its neighbours. Gradient vector flow diffuses the edge-map gradient into the flat background, so its force is still about half its maximum five pixels away after 128 iterations, and ten pixels away after 1,000, when more than a third of it remains twenty pixels away.

The two snakes with the traditional edge force started outside both reach the object, because tension alone shrinks a closed snake until it meets the edges (Cohen’s rectangle example, below, describes the same behaviour). Both bridge the notch: they leave out 6% and 16% of it. The NumPy snake’s IoU of 0.866 is lost almost entirely in the notch. The same force started from a small circle inside the object, where the smoothed image is nearly flat, has nothing to stop its contraction and collapses to the minimum of 8 points, with an IoU of 0; its notch score of 1 is meaningless here and is why the two columns must be read together.

scikit-image’s default run is the clearest lesson about stopping rules. It ends almost where it started (IoU 0.452, close to the area ratio of the initial circle), because its convergence test stops the loop as soon as the snake lies within 0.1 pixels of any one of its last ten saved positions, the most recent of which is one iteration old, and a snake creeping inward under weak tension moves less than that in one iteration (our reading of the code). With convergence=0 it runs all 2,500 iterations and behaves like the NumPy snake: IoU 0.833, the notch bridged. A snake that stops is not necessarily a snake that has converged.

Balloon forces depend on how the image force is normalized, and the two sources normalize differently. Cohen (1991) normalizes the edge force to a fixed magnitude $k$ in the direction of $-\nabla P$, so that weak and strong gradients pull alike, and chooses $k$ slightly larger than the pressure $k_1$; Xu and Prince (1998) scale the external forces in all their experiments to a maximum magnitude of one. In this noisy image the smoothed noise has a gradient almost everywhere, and Cohen’s normalization gives every noise gradient the full magnitude $k = 0.5$, larger than the pressure. With it, neither inflating balloon matches the object (IoU 0.354 and 0.183), and the two deflating balloons leave the notch out (fractions 0.990 and 0.974) but score lower IoU than the snakes that bridge it. The scores do not show where these balloons stopped; our reading is that noise gradients at full strength can stop a weak balloon, and that a pressure close to $k$ can also carry it across the true edges. The same pressure added to the edge force scaled by its maximum, as Xu and Prince scale their forces, gives the behaviour Cohen describes: from inside or outside the balloon runs to the boundary and stops there, with IoU 0.985 and 0.971 and the notch respected. The inflating balloon respects the notch naturally, since seen from inside the notch is a bump in the boundary, not a cavity. Which normalization suits an image depends, in our reading, on its noise; the example shows both because the sources differ.

Gradient vector flow behaves as Xu and Prince report. After 1,000 diffusion iterations the snake follows the notch whether it starts outside or inside (IoU 0.982, about 97% of the notch left out), without being told which way to move. After 128 iterations, the $\sqrt{N}$ that Xu and Prince quote in their timing, applied to this image, the field has not yet spread into the notch strongly enough: the snake enters part of the way (37% of the notch left out). How far the field reaches depends on how long the diffusion runs; Xu and Prince define the field as the steady state of the diffusion and give an iteration count only with their timing.

Two cases were added after the first run and are disclosed here. The scikit-image run without early stopping was added after the default run was seen to stop within a few pixels of its start. The two balloons on the maximum-scaled edge force were added after the balloons on Cohen’s normalized force scored as printed above. Two repairs that did not change what the code computes were also made before the final run: a small constant in the denominator of the normal, after two coincident points produced a division by zero, and a sparse LU factorization in place of a dense inverse, after the growing point count of one balloon run made the dense inverse slow.

Properties and Behavior

Local minima by design. The energy is not convex, and Kass et al. did not want a global minimum: their wood-grain example (Fig. 1) shows three different local minima, three plausible organizations of two dark lines, and they argue that choosing among such alternatives can need knowledge that low-level processing lacks. The snake finds the minimum nearest its start, so the start matters.

Capture range. The distance from which an edge can attract the snake is set by the reach of the external force. For the traditional potential this is a few times $\sigma$. Kass et al. enlarged it by scale-space continuation: let the snake reach equilibrium on a heavily blurred potential, which attracts it from far away but localizes the edge poorly, then reduce the blur step by step (their Fig. 4). Xu and Prince (1998) list multiresolution methods, pressure forces and distance potentials as the main earlier ways of enlarging the capture range.

Shrinking. With no image force the internal energy has its minimum at a point (closed snake) or a straight line (open snake with fixed ends), so a snake started in a flat region contracts. Cohen (1991) illustrated this with a curve inside a black rectangle on a white background: even with perfect edge detection, the curve vanishes unless some of its points start close enough to an edge to be caught; a curve surrounding the rectangle, by contrast, shrinks onto it. Isolated noise maxima of the gradient can also stop a passing snake.

Concavities. Xu and Prince (1998, Fig. 1) showed that inside a U-shaped object’s concavity the traditional forces point horizontally in opposite directions, toward the two walls, so the snake is pulled apart instead of drawn down, and wrote that no choice of $\alpha$ and $\beta$ corrects this. Distance potential forces (Cohen & Cohen, 1993, as Xu and Prince describe them) enlarge the capture range but point the same way in the concavity and fail likewise. Xu and Prince believe the failure is a property of all conservative forces, those that are the gradient of a potential.

Smoothing of corners. The internal energy rounds sharp corners unless $\beta$ is set to zero there; Cohen’s drawn test figure (his Fig. 3) and Xu and Prince’s room-shaped test figure (their Fig. 5) show slightly rounded corners, which both attribute to the regularization.

Hysteresis and tracking. Because a snake keeps minimizing as the image changes, it follows a slowly moving feature and lets go only when the feature moves too far, then reattaches when the feature comes back (Kass et al., Fig. 6). Kass et al. tracked a speaker’s lips through a two-second video with edge snakes initialized by hand in the first frame and no constraints between frames (Fig. 8), and suggested giving the snake mass so it would predict its next position from its velocity. They also coupled left and right snakes with a disparity-smoothness term for stereo matching (Eq. 9). Later contour trackers, such as CONDENSATION, also follow outlines represented as parametric curves, there B-splines.

Not region-aware. A classical snake sees the image only along its own path, through edge, line and termination terms. Region statistics inside and outside the curve enter only in later models; McInerney and Terzopoulos (1996) list several groups that added region-based information to make snakes less sensitive to spurious edges and to the starting position.

Limitations

McInerney and Terzopoulos (1996) summarize the limits of snakes as they stood in the mid-1990s:

• Designed for interaction. Snakes were built as interactive tools; used without a user, they must be initialized close to the structure to work well.
• Geometric inflexibility. The internal energy keeps the snake from representing long thin tubes, protrusions or branches.
• Fixed topology. A parametric snake is one curve. It cannot split to outline two objects or merge two curves, so the topology of the structure must be known in advance. Caselles, Kimmel and Sapiro (1997) add that a snake started around several objects ends, without special procedures, near a figure like their convex hull.
• Parameterization dependence. The energy changes with the speed at which the curve is traversed, so it is not a property of the curve’s shape alone; Caselles et al. call it non-intrinsic.

The example adds practical limits. Results depend on several coupled weights, on the force normalization and on the stopping rule; a snake that stops early can look converged. Points bunch up or spread apart as the curve moves, so implementations resample, and a resampling rule is one more choice the original paper does not make. Balloon forces need a direction chosen in advance, and pressure strong enough to cross noise can also cross the true boundary (Xu & Prince, 1998, citing Tek and Kimia on weak edges being overwhelmed).

Against later segmentation methods, the snake optimizes locally from an initial guess, while graph cuts find a global optimum of a two-label energy of the kind they can represent and the watershed partitions the whole image without a starting curve; learned segmenters such as U-Net replace hand-designed forces altogether. Active contours remain a way to refine a boundary with sub-pixel precision and a smoothness prior when a good starting curve is available.

Variants

Balloons

Cohen (1991) addressed both the instability of the explicit image step and the shrinking. He first replaced the force $-\nabla P$ by a force of fixed magnitude in the same direction,

$$\mathbf{F} = -k \frac{\nabla P}{\lVert \nabla P \rVert},$$

with $k$ about the size of a pixel, so that no step jumps across an edge. He then made the curve behave like a balloon being inflated by adding a pressure along the unit normal $\mathbf{n}(s)$:

$$\mathbf{F} = k_1 \mathbf{n}(s) - k \frac{\nabla P}{\lVert \nabla P \rVert} .$$

Changing the sign of $k_1$, or the orientation of the curve, deflates instead. Cohen chose $k$ and $k_1$ of the same order, below a pixel, with $k$ slightly larger so that an edge can stop the inflation. A weak edge, or an isolated noise point, makes a kink in the expanding curve that smoothing and pressure remove after a few iterations, so the balloon passes it. At equilibrium the curve sits slightly outside the true contour, because the edge force balances the pressure and the smoothing; reducing the pressure then localizes it. Cohen also built potentials from a binary edge map, from the Canny–Deriche detector, convolved with a Gaussian, so that a snake could close a broken set of Canny edges, and he applied balloons to ultrasound and MR images of the heart and to slice-by-slice reconstruction of the ventricles in 3D.

Xu and Prince (1998) note the costs: a balloon must be told whether to push out or in, may overwhelm weak boundaries if it pushes too hard, and may not enter concavities if it pushes the wrong way or too weakly. In their room-shaped test figure (Fig. 5b), a balloon strong enough to reach the alcoves also bulged out through the open doorways that a snake should have bridged. Cohen and Cohen (1993) developed balloons with finite elements for 2D and 3D images; McInerney and Terzopoulos (1996) note that inflation forces also appear in Terzopoulos et al. (1988).

Gradient vector flow

Xu and Prince (1998) replaced the potential force by a vector field $\mathbf{v}(x, y) = (u(x, y), v(x, y))$ that minimizes (their Eq. 12)

$$\mathcal{E} = \iint \mu \left( u_x^2 + u_y^2 + v_x^2 + v_y^2 \right) + \lvert \nabla f \rvert^2 \lvert \mathbf{v} - \nabla f \rvert^2 \, dx\, dy ,$$

where $f$ is an edge map, larger near edges, for example $\lvert \nabla (G_\sigma * I) \rvert^2$, normalized to $[0, 1]$. Near edges, where $\lvert \nabla f \rvert$ is large, the second term keeps $\mathbf{v}$ close to $\nabla f$; in flat regions the first term makes the field vary slowly. The smoothness term is the one Horn and Schunck used for optical flow. The minimizer satisfies

$$\mu \nabla^2 u - (u - f_x)(f_x^2 + f_y^2) = 0, \qquad \mu \nabla^2 v - (v - f_y)(f_x^2 + f_y^2) = 0 ,$$

two decoupled linear equations (Eq. 13). In a homogeneous region the second term vanishes, so $u$ and $v$ obey Laplace’s equation and are interpolated from the surrounding boundary, which Xu and Prince call a competition among boundary vectors and which yields vectors pointing into concavities. They solve the equations as diffusion to steady state with an explicit scheme (Eq. 16), stable when $\mu \Delta t / (\Delta x \Delta y) \le 1/4$ (Eq. 18), and set $\mu$ according to the noise (more noise, larger $\mu$).

The GVF snake is $\mathbf{x}_t = \alpha \mathbf{x}'' - \beta \mathbf{x}'''' + \mathbf{v}$ (Eq. 10). Because $\mathbf{v}$ is not in general the gradient of a potential, the snake no longer minimizes an energy; it is defined by force balance. Xu and Prince showed on line drawings, a noisy U shape and an MR image of the left ventricle that GVF snakes can start inside, outside or across the boundary, have a large capture range, and enter concavities, and they extended GVF to 3D deformable surfaces. They left open a complete characterization of the capture range and how to choose $\mu$.

Other optimizers and models

The Euler–Lagrange descent is not the only way to minimize a snake energy. McInerney and Terzopoulos (1996) list dynamic programming (Amini et al., 1990), for a more extensive search toward global minima, and simulated annealing (Poon et al., 1994; Grzeszczuk and Levin, 1994). Topologically adaptable snakes, or T-snakes (McInerney and Terzopoulos, 1995, as the survey reports them), can sense and change their topology: in the survey’s vertebra example a snake starts as one closed curve and ends as three. Deformable surfaces extend the same ideas to 3D volumes (McInerney & Terzopoulos, 1996).

Geometric and implicit contours

A different line of work dropped the parameterized curve. Caselles, Catté, Coll and Dibos (1993) and Malladi, Sethian and Vemuri (1994, 1995), as Caselles, Kimmel and Sapiro (1997) describe them, moved the curve by a geometric flow, a speed along the normal made of a constant term and a curvature term, scaled by an edge-stopping function $g$ such as $1/(1 + \lvert \nabla \hat{I} \rvert^p)$ for a smoothed image $\hat{I}$. Represented as the zero level set of a function on the image grid, following Osher and Sethian (1988), the curve changes topology automatically and several objects can be detected without knowing how many there are. Caselles, Kimmel and Sapiro note that the constant speed term is similar to Cohen’s balloon force.

Geodesic active contours (Caselles, Kimmel & Sapiro, 1997) connected the two lines. Starting from a snake with $\beta = 0$, and generalizing the edge term to $g(\lvert \nabla I \rvert)$, they used Maupertuis’s principle to show that minimizing the snake energy, with a particular choice of a free energy constant, is equivalent to finding a curve of minimal weighted length

$$L_R = \int_0^{L(C)} g\left( \lvert \nabla I(C(s)) \rvert \right) ds ,$$

a geodesic in a metric defined by the image (their Eqs. 8 and 12). The weighted length does not depend on how the curve is parameterized. Its gradient descent moves each point along the inward normal $\mathbf{N}$ with speed $g\kappa - \nabla g \cdot \mathbf{N}$ (Eq. 13), where $\kappa$ is the curvature: curvature smoothing weighted by the edge function, plus a term that pulls the curve toward the bottom of the edge valley. The second term is what earlier geometric models lacked; the authors report that it lets the contour follow boundaries whose gradient varies strongly, including small gaps. The curvature flow replaces the snake’s bending term as the regularizer, which is why they could set $\beta = 0$.

Implementing these flows on a grid, with upwind schemes, reinitialization and narrow bands, along with region-based models such as Chan–Vese, is the subject of Level Set Methods. scikit-image 0.25 provides chan_vese and morphological approximations, morphological_geodesic_active_contour and morphological_chan_vese, beside the parametric active_contour. Xu and Prince suggested that GVF might also serve as an external force for geometric snakes.