Superpixels

Superpixels are small, connected regions of similar pixels that replace the pixel grid as the unit of later processing. A superpixel algorithm deliberately oversegments an image into hundreds or thousands of pieces, each meant to lie inside one object, so that later steps work on hundreds of regions instead of hundreds of thousands of pixels. Superpixels are a preprocessing tool for image segmentation and other tasks.

Definition

A superpixel segmentation is an oversegmentation: a partition of the image into many connected regions that are uniform in color or texture and adhere to image boundaries, so that almost every true object boundary is also a superpixel boundary. Superpixels are usually also expected to be compact, of similar size, fast to compute and controllable in number; Stutz, Hermans and Leibe (2018) use the last requirement to separate superpixel algorithms from segmentation algorithms run in “oversegmentation mode”.

Ren and Malik (2003) introduced the term. They grouped pixels into superpixels with normalized cuts before learning which groupings form good segments, arguing that pixels are not natural entities, only artifacts of discretization, and too numerous for pixel-level optimization. Measuring how many human-marked boundaries survived, they found 200 superpixels sufficient for their 240 × 160 images.

Intuition

Neighboring pixels inside an object are largely redundant, so merging them early loses little. The bet is that no superpixel straddles a boundary that matters: no later step can recover a boundary inside one, so the algorithm errs toward too many regions.

A competing wish is regularity. Like the pixel grid they replace, compact cells of similar size have few neighbors and comparable features. Boundary adherence pulls cells into irregular shapes; regularity pulls them back toward a grid. Many algorithms expose this trade-off as a parameter.

Formal Definition

Given an image on a pixel set $\Omega$ of size $N$ , a superpixel segmentation is a partition $\mathcal{S} = \{S_1, \dots, S_K\}$ of $\Omega$ into connected sets, with $K \ll N$ . It is judged against a ground-truth segmentation $\mathcal{G} = \{G_1, \dots, G_M\}$ drawn by people, as in the Berkeley Segmentation Dataset, mainly with two measures.

Boundary recall is the fraction of ground-truth boundary pixels that lie close to some superpixel boundary:

\mathrm{Rec}(\mathcal{G}, \mathcal{S}) = \frac{TP}{TP + FN},

where a ground-truth boundary pixel is a true positive if a superpixel boundary lies within a small tolerance of it. Stutz et al. credit the measure to Martin, Fowlkes and Malik (2004), who used it to evaluate boundary detectors. Achanta et al. (2012) used a tolerance of two pixels; Stutz et al. used a $(2r+1) \times (2r+1)$ neighborhood with $r$ equal to 0.0025 times the image diagonal, rounded, which gives $r = 1$ on BSDS500.

Undersegmentation error measures how far superpixels leak across ground-truth regions. Stutz et al. trace it to Levinshtein et al. (2009), whose form charges every superpixel touching a ground-truth region in full, even if it overlaps by one pixel; Achanta et al. counted only superpixels with more than 5% of their pixels inside the region. Stutz et al. adopted a form they attribute to Neubert and Protzel, which charges each overlapping superpixel the smaller of its parts inside and outside the region:

\mathrm{UE}(\mathcal{G}, \mathcal{S}) = \frac{1}{N} \sum_{i=1}^{M} \sum_{S_j \cap G_i \neq \varnothing} \min\big(|S_j \cap G_i|,\; |S_j \setminus G_i|\big).

Regularity is scored separately, for example by the compactness of Schick et al., as given by Stutz et al.: the area-weighted average of $4\pi A(S_j) / P(S_j)^2$ , which is 1 for a continuous disk of area $A$ and perimeter $P$ .

Families of Algorithms

Achanta et al. divided superpixel methods into graph-based and gradient-ascent methods; Stutz et al. used finer categories for the 28 algorithms they evaluated.

Graph-based. Pixels are nodes and edges carry similarities. Normalized cuts superpixels, as used by Ren and Malik, recursively partition this graph; Achanta et al. found them very regular but poor at boundary adherence, and the slowest method they compared. Felzenszwalb and Huttenlocher (2004) process edges in order of increasing weight and join two components unless the edge between them exceeds either component’s internal difference plus a threshold $k/|C|$ . The method runs in $O(N \log N)$ and adheres well to boundaries, but gives very irregular regions and no direct control of their number. Entropy rate superpixels (Liu et al., 2011) select edges greedily to maximize the entropy rate of a random walk on the graph, which favors compact, homogeneous clusters, plus a term that balances sizes.

Gradient ascent and mode seeking. Watershed (Vincent & Soille, 1991) floods the gradient image from its minima; mean shift (Comaniciu & Meer, 2002) and quick shift (Vedaldi & Soatto, 2008) assign each pixel to a mode of a density over color and position. None directly controls the number or compactness of regions. TurboPixels (Levinshtein et al., 2009) dilate regularly spaced seeds by a level-set geometric flow guided by image gradients, giving compact superpixels of uniform size.

Clustering. SLIC (Achanta et al., 2012) runs a localized k-means on color and position, as described below.

Energy optimization. SEEDS (Van den Bergh et al., 2015) starts from a regular grid and refines it by hill climbing: it moves blocks of pixels, then single boundary pixels, between neighboring superpixels when the move improves an energy that rewards concentrated color histograms and smooth boundaries.

Learned superpixels. Superpixel Sampling Networks (Jampani et al., 2018) make SLIC differentiable: hard nearest-center assignments become soft associations $Q_{pi} = e^{-\lVert I_p - S_i \rVert^2}$ with the surrounding centers, and the clustered features come from a convolutional network, concatenated with position and Lab color. Trained end to end with a task-specific reconstruction loss and a compactness loss, the network can learn superpixels suited to semantic segmentation or optical flow.

Stutz et al. optimized each algorithm’s parameters on training sets, enforced connectivity, and ranked algorithms by boundary miss rate ( $1 - \mathrm{Rec}$ ) plus undersegmentation error, each averaged over 200 to 5,200 superpixels, on five datasets. They recommended six for practice: ETPS, SEEDS, ERS, CRS, ERGC and SLIC.

SLIC in Detail

SLIC (simple linear iterative clustering) clusters pixels in a five-dimensional space of CIELAB color $[l, a, b]$ and position $[x, y]$ ; CIELAB is the color space in which Euclidean distance roughly tracks perceived color difference. As described by Achanta et al. (2012):

Initialization. For $N$ pixels and $K$ desired superpixels, place $K$ centers on a regular grid with step $S = \sqrt{N/K}$ , then move each to the lowest-gradient position in its $3 \times 3$ neighborhood, so that no center sits on an edge or a noisy pixel.
Assignment. Each center compares itself only with pixels in a $2S \times 2S$ window around it, and each pixel takes the label of the closest center that reaches it.
Update. Each center moves to the mean $[l, a, b, x, y]$ of its pixels. The authors found 10 iterations of assignment and update enough for most images.
Connectivity. SLIC does not guarantee connected clusters, so a final pass reassigns stray pieces to nearby superpixels.

The distance between pixel $i$ and center $k$ combines color and spatial distances,

d_c = \lVert [l, a, b]_i - [l, a, b]_k \rVert, \qquad d_s = \lVert [x, y]_i - [x, y]_k \rVert, \qquad D = \sqrt{d_c^2 + \left(\frac{d_s}{S}\right)^2 m^2}.

Dividing $d_s$ by the grid step makes the balance independent of superpixel size, and the compactness $m$ stands in for the largest color distance expected within a cluster. Large $m$ gives compact superpixels; small $m$ gives tighter boundary adherence with less regular shapes. The authors suggest $m$ in $[1, 40]$ for CIELAB, with 10 as the default. Weighing color against position resembles the range and domain kernels of the bilateral filter.

A pixel falls within the windows of fewer than eight centers, so SLIC is $O(N)$ , independent of $K$ , whereas k-means over the whole image costs $O(NKI)$ for $I$ iterations. Achanta et al. reported SLIC as the fastest method they compared ( $2048 \times 1536$ pixels in 14.94 s on a 2.26 GHz machine), with the lowest undersegmentation error on the Berkeley dataset over most superpixel counts, and second to Felzenszwalb–Huttenlocher in boundary recall, or first with a smaller $m$ .

Examples

SLIC on a synthetic image with known regions

The example builds a $96 \times 96$ sRGB image with four ground-truth regions (a background with a left-to-right color ramp, a red disk, a blue square, and a pale triangle only 4 to 17 CIELAB units from the background beside it) and adds Gaussian noise. It runs SLIC in NumPy and scores it with boundary recall at tolerances of one and two pixels, undersegmentation error and compactness, against a grid of squares with the same number of cells. Three details are this implementation’s choices: the seed-moving gradient is the squared central-difference gradient of the Lab image, each stray piece goes to the neighbor it shares the longest border with, and perimeters are counted in pixel edges.

import numpy as np
from scipy import ndimage

# A 96 x 96 test image with known regions: a background whose color changes
# smoothly from left to right, a red disk, a blue square, and a pale triangle
# whose color differs only slightly from the background around it.
H = W = 96
yy, xx = np.mgrid[0:H, 0:W]
gt = np.zeros((H, W), int)                                   # ground-truth region ids
gt[(yy - 30)**2 + (xx - 32)**2 <= 18**2] = 1                 # disk
gt[52:84, 50:84] = 2                                         # square
gt[(yy >= 10) & (yy <= 44) & (np.abs(xx - 72) <= (yy - 10) / 2)] = 3   # triangle
colors = np.array([[0, 0, 0], [0.80, 0.20, 0.20], [0.20, 0.30, 0.80], [0.62, 0.70, 0.68]])
t = xx / (W - 1)
background = np.stack([0.85 - 0.35 * t, 0.80 - 0.10 * t, 0.55 + 0.25 * t], axis=-1)
rgb = np.where((gt == 0)[..., None], background, colors[gt])


def rgb_to_lab(rgb):
    """sRGB in [0, 1] -> CIELAB (D65 white), as defined by IEC 61966-2-1 and CIE 015."""
    c = np.where(rgb <= 0.04045, rgb / 12.92, ((rgb + 0.055) / 1.055)**2.4)
    M = np.array([[0.4124, 0.3576, 0.1805],
                  [0.2126, 0.7152, 0.0722],
                  [0.0193, 0.1192, 0.9505]])
    xyz = c @ M.T / np.array([0.95047, 1.0, 1.08883])
    f = np.where(xyz > (6 / 29)**3, np.cbrt(xyz), xyz / (3 * (6 / 29)**2) + 4 / 29)
    return np.stack([116 * f[..., 1] - 16,
                     500 * (f[..., 0] - f[..., 1]),
                     200 * (f[..., 1] - f[..., 2])], axis=-1)


def slic(lab, K, m=10.0, iters=10):
    h, w, _ = lab.shape
    S = np.sqrt(h * w / K)                                   # grid interval
    # 1. Centers on a regular grid, moved to the lowest-gradient pixel in 3 x 3.
    gy, gx = np.gradient(lab, axis=(0, 1))
    grad = (gx**2 + gy**2).sum(-1)
    ys = np.arange(S / 2, h, S).astype(int)
    xs = np.arange(S / 2, w, S).astype(int)
    centers = []
    for y in ys:
        for x in xs:
            win = grad[max(y - 1, 0):y + 2, max(x - 1, 0):x + 2]
            dy, dx = np.unravel_index(np.argmin(win), win.shape)
            y2, x2 = max(y - 1, 0) + dy, max(x - 1, 0) + dx
            centers.append([*lab[y2, x2], y2, x2])
    C = np.array(centers, float)                             # rows: l, a, b, y, x
    R = int(np.ceil(S))                                      # 2S x 2S search window
    for _ in range(iters):
        # 2. Assignment: each center competes only for pixels near it.
        dist = np.full((h, w), np.inf)
        label = np.full((h, w), -1)
        for k, (l, a, b, cy, cx) in enumerate(C):
            y0, y1 = max(int(cy) - R, 0), min(int(cy) + R + 1, h)
            x0, x1 = max(int(cx) - R, 0), min(int(cx) + R + 1, w)
            dc2 = ((lab[y0:y1, x0:x1] - (l, a, b))**2).sum(-1)
            ds2 = (yy[y0:y1, x0:x1] - cy)**2 + (xx[y0:y1, x0:x1] - cx)**2
            D = np.sqrt(dc2 + ds2 / S**2 * m**2)
            better = D < dist[y0:y1, x0:x1]
            dist[y0:y1, x0:x1][better] = D[better]
            label[y0:y1, x0:x1][better] = k
        # 3. Update: each center moves to the mean labxy of its pixels.
        feats = np.concatenate([lab, yy[..., None], xx[..., None]], -1).reshape(-1, 5)
        n = np.bincount(label.ravel(), minlength=len(C))
        for j in range(5):
            s = np.bincount(label.ravel(), feats[:, j], minlength=len(C))
            C[n > 0, j] = s[n > 0] / n[n > 0]
    return enforce_connectivity(label)


def enforce_connectivity(label):
    """Keep the largest connected piece of each superpixel; give every other
    piece to the neighboring superpixel it shares the longest border with."""
    out = label.copy()
    orphans = []
    for k in np.unique(label):
        comp, n = ndimage.label(label == k)
        if n > 1:
            sizes = np.bincount(comp.ravel())[1:]
            orphans += [comp == i + 1 for i in range(n) if i != np.argmax(sizes)]
    for piece in orphans:
        ring = ndimage.binary_dilation(piece) & ~piece
        out[piece] = np.bincount(out[ring]).argmax()
    _, out = np.unique(out, return_inverse=True)             # relabel 0..n-1
    return out.reshape(label.shape)


def boundaries(seg):
    b = np.zeros(seg.shape, bool)
    b[:, :-1] |= seg[:, :-1] != seg[:, 1:]
    b[:-1, :] |= seg[:-1, :] != seg[1:, :]
    return b


def boundary_recall(seg, gt, r=2):
    """Fraction of ground-truth boundary pixels within r pixels of a superpixel boundary."""
    near = ndimage.binary_dilation(boundaries(seg), np.ones((2 * r + 1, 2 * r + 1)))
    g = boundaries(gt)
    return (near & g).sum() / g.sum()


def undersegmentation_error(seg, gt):
    """Neubert and Protzel's form: each superpixel overlapping a ground-truth
    region pays the smaller of its part inside and its part outside that region."""
    overlap = np.zeros((seg.max() + 1, gt.max() + 1))
    np.add.at(overlap, (seg.ravel(), gt.ravel()), 1)
    size = overlap.sum(1, keepdims=True)
    return np.minimum(overlap, size - overlap)[overlap > 0].sum() / seg.size


def compactness(seg):
    """Schick et al.'s compactness: area-weighted 4*pi*A/P^2, perimeter counted in pixel edges."""
    p = np.pad(seg, 1, constant_values=-1)
    total = 0.0
    for k in range(seg.max() + 1):
        mask = p == k
        P = (mask[:, 1:] != mask[:, :-1]).sum() + (mask[1:, :] != mask[:-1, :]).sum()
        A = mask.sum()
        total += A * 4 * np.pi * A / P**2
    return total / seg.size


def grid(n_target):
    side = int(round(np.sqrt(n_target)))                     # side x side squares
    return (yy * side // H) * side + (xx * side // W)


def report(name, seg, gt):
    print(f"{name:22s} n={seg.max() + 1:4d}  recall(r=1) {boundary_recall(seg, gt, 1):.3f}  "
          f"recall(r=2) {boundary_recall(seg, gt, 2):.3f}  UE {undersegmentation_error(seg, gt):.3f}  "
          f"CO {compactness(seg):.3f}")


rng = np.random.default_rng(0)
img = np.clip(rgb + rng.normal(0, 0.03, rgb.shape), 0, 1)    # noisy sRGB image
lab = rgb_to_lab(img)
for K in (50, 150):
    report(f"SLIC K={K}, m=10", slic(lab, K, m=10), gt)
    report(f"grid {int(round(np.sqrt(K)))**2} squares", grid(K), gt)
for m in (1, 5, 20, 40):
    report(f"SLIC K=50, m={m}", slic(lab, 50, m=m), gt)

Output:

SLIC K=50, m=10        n=  49  recall(r=1) 0.944  recall(r=2) 0.971  UE 0.044  CO 0.326
grid 49 squares        n=  49  recall(r=1) 0.471  recall(r=2) 0.627  UE 0.214  CO 0.785
SLIC K=150, m=10       n= 144  recall(r=1) 0.995  recall(r=2) 1.000  UE 0.017  CO 0.366
grid 144 squares       n= 144  recall(r=1) 0.532  recall(r=2) 0.802  UE 0.181  CO 0.785
SLIC K=50, m=1         n=  49  recall(r=1) 0.979  recall(r=2) 1.000  UE 0.087  CO 0.146
SLIC K=50, m=5         n=  49  recall(r=1) 0.987  recall(r=2) 1.000  UE 0.049  CO 0.242
SLIC K=50, m=20        n=  49  recall(r=1) 0.923  recall(r=2) 0.958  UE 0.052  CO 0.514
SLIC K=50, m=40        n=  49  recall(r=1) 0.812  recall(r=2) 0.876  UE 0.087  CO 0.619

Asked for 50 and 150 superpixels, SLIC returned 49 and 144, from seed grids of $7 \times 7$ and $12 \times 12$ . With $m = 10$ and 49 superpixels, it recovered 94.4% of ground-truth boundary pixels within one pixel, against 47.1% for 49 squares, with undersegmentation error 0.044 against 0.214. The squares score 0.785 in compactness, $\pi/4$ , the most any shape can score when perimeters are counted in pixel edges, but they cut through every object. With 144 cells, the grid’s two-pixel recall reaches 0.802 because dense boundaries land near true ones by chance.

The sweep over $m$ shows the trade-off. Compactness rises from 0.146 at $m = 1$ to 0.619 at $m = 40$ , while one-pixel recall falls from about 0.98 to 0.812, all of the loss on the triangle’s weak boundary. Undersegmentation error is lowest for $m$ from 5 to 20 (0.044 to 0.052) and reaches 0.087 at both ends, mostly from leakage between triangle and background. At $m = 40$ compact superpixels ignore the weak edge; at $m = 1$ , in our reading, color dominates the distance and noise makes outlines ragged (60% more boundary pixels than at $m = 10$ ). Over noise seeds 0 to 9, $m = 10$ gave one-pixel recall 0.899 to 0.960 and undersegmentation error 0.031 to 0.069; the compactness and recall trends held for every seed, and the error was lowest at an intermediate $m$ for 9 of the 10.

Common Misconceptions

“Superpixels are a segmentation.” They know nothing of objects or classes; objects are recovered by later steps that merge or label superpixels.
“Higher boundary recall means better superpixels.” Recall rises with more superpixels or more ragged boundaries, so it is compared at equal counts and read with undersegmentation error.
“Superpixels are always connected.” Clustering methods, SLIC included, need a connectivity pass. Achanta et al. noted that the quick shift code they used did not guarantee connected superpixels, and Stutz et al. enforced connectivity for every algorithm.

Where It Is Used

Achanta et al. used superpixels as the units of multiclass segmentation on MSRC and PASCAL VOC 2010, and to make graph-based segmentation of large electron-microscopy images and volumes tractable. Selective search (Uijlings et al., 2013), the proposal method of the R-CNN detector, starts from Felzenszwalb–Huttenlocher regions and merges them hierarchically into object proposals. Stutz et al. list uses in tracking, stereo, depth and optical-flow estimation, saliency and dataset annotation. Watershed superpixels come from mathematical morphology.