Shape Descriptors

A shape descriptor is a compact description of the outline or silhouette of an object, computed after the object has been separated from its background, that can be compared between objects regardless of where they are, how large they are and, usually, how they are turned. Shape descriptors range from single numbers measured on a region, such as its area, circularity or solidity, through codes and series that encode the boundary curve, such as Freeman chain codes and Fourier descriptors, to point-wise histograms that support matching between deformed shapes, such as the shape context of Belongie, Malik and Puzicha (2002). This article defines these families, shows how each is normalized for translation, scale, rotation and the starting point on the contour, and measures in code what the normalizations achieve and what they lose. Moment-based descriptors, which summarize a region by weighted sums of its pixel coordinates (Hu’s invariants, Zernike moments), are covered in image moments; this article refers to them only in summary.

Shape was one of the global cues, with colour and texture, by which early image retrieval systems described whole images. In our reading it matters most where colour and texture carry little information, as in the silhouettes, handwritten digits and trademarks of Belongie et al.’s experiments.

Definition

A shape descriptor maps a shape, given as a binary region, a closed boundary curve or a set of points sampled from the object’s edges, to a vector or a set of vectors, together with a way of comparing two descriptors. It is useful when it is unchanged by the transformations that do not change the shape (at least translation and uniform scaling, often rotation, sometimes reflection) and changes when the shape itself changes.

Three families are distinguished by what they start from:

• Region descriptors are computed from all pixels of a region: area, perimeter, circularity, eccentricity, solidity, extent, the Euler number, and image moments.
• Contour descriptors are computed from the ordered boundary curve: chain codes, curvature and other signatures, and Fourier descriptors.
• Point-set descriptors are computed from unordered points on internal and external edges, such as the shape context, and compare shapes through correspondences between points.

A further distinction is between descriptors that give one vector for the whole shape and those that give one vector per boundary point or part, which are then matched point by point.

Intuition

Imagine describing a cut-out shape to someone over the telephone. A few numbers give the gist: how big it is, how round, whether it has dents or holes. A precise description walks around the outline: “go right three steps, up one, diagonally up two”, which is a chain code. A smoother description says the outline is a circle, plus a slight oval, plus a three-lobed ripple, plus finer ripples, which is a Fourier description; dropping the finest ripples gives a recognizable but rounded version of the shape. A third way describes the view from each point of the outline: where the rest of the outline lies, near or far, in which direction. If two shapes are similar, each point on one has a point on the other with a similar view, and pairing up the points shows how one shape deforms into the other; this is the idea of the shape context.

Each description has to forget something. The walking directions depend on where the walk starts and which way the outline faces; the Fourier description depends on the same, but in a way that can be removed exactly; the view from each point depends on which direction is called “up”. Each normalization removes one of these dependences, and each removal can also erase a distinction that mattered: a 6 and a 9, or a shape and its mirror image.

Formal Definition

Region descriptors

For a binary region with area $A$ and perimeter $P$, common scalar descriptors are listed below, with their definitions in scikit-image’s regionprops (version 0.25.2), which computes all of them:

• Circularity, also called compactness, $4\pi A / P^2$, which is 1 for a disk and smaller for every other shape. scikit-image does not provide it directly; it is formed from area and perimeter.
• Eccentricity of the ellipse with the same second moments as the region, from 0 for a circle towards 1 for a line; it and the orientation and axis lengths are second-order moments, derived in the Image Moments article.
• Solidity, the ratio of the region’s pixels to the pixels of its convex hull, which is below 1 when the region has concavities.
• Extent, the ratio of the region’s pixels to those of its bounding box.
• Euler number, the number of connected components minus the number of holes, a topological count that does not change under any continuous deformation.
• Maximum Feret diameter, the longest distance between points of the region’s convex hull.

Area, solidity and the Euler number are counts of pixels and are well defined on a grid. The perimeter is not: a digital region has no unique boundary length, and each estimator approximates it differently. scikit-image’s perimeter approximates the contour, in its own words, as “a line through the centers of border pixels”, weighting the steps between them, and perimeter_crofton uses the Crofton formula of integral geometry, which counts how many times lines in 2 or 4 directions cross the boundary. The first example below measures how both, and the chain-code length, depend on the region’s orientation. Any error in $P$ enters circularity squared.

Mathematical morphology supplies a further family of region descriptors through the skeleton, or medial axis, which reduces a shape to thin curves along its middle that preserve its topology.

Chain codes

A boundary traced on the pixel grid is a closed path whose consecutive pixels are neighbours. The Freeman chain code (OpenCV’s source calls its reader for such codes a “Freeman chain reader”) records the direction of each step as one of eight symbols: with $x$ to the right and $y$ up, code $k \in \{0, \dots, 7\}$ is the step at angle $45° \cdot k$, so 0 is east, 2 north, 4 west and 6 south, and odd codes are diagonal. The code, together with the start pixel, reproduces the boundary exactly. OpenCV’s legacy C interface defines a chain-code mode (CV_CHAIN_CODE) and a Freeman chain reader, and its contour finder can return every boundary pixel (CHAIN_APPROX_NONE) or only the end points of horizontal, vertical and diagonal runs (CHAIN_APPROX_SIMPLE).

The code depends on the start pixel and on the orientation. Two normalizations follow from its definition:

• Rotation by multiples of 90°. Rotating the image by a quarter turn adds 2 (mod 8) to every symbol; the first difference $d_i = (c_i - c_{i-1}) \bmod 8$, which records the turn at each step, is unchanged. Adding 1 to every symbol would be the analogue for 45°, but it does not describe a rotated shape: it turns axis steps of length 1 into diagonal steps of length $\sqrt 2$, and a shape rasterized at 45° is a different set of pixels, not a relabelled one.
• Start point. The difference code is a cyclic sequence; choosing the cyclic rotation that is smallest in lexicographic order gives a representation that does not depend on the start pixel. The code below calls it the shape number.

A chain code also gives a length estimate: each even step has length 1 and each odd step $\sqrt 2$, so the boundary length is $n_{\text{even}} + \sqrt 2\, n_{\text{odd}}$. For a digital straight edge at angle $\alpha$ between 0° and 45° to the grid, the steps are a mix of axis and diagonal moves whose total length per unit length of edge is $\cos\alpha + (\sqrt 2 - 1)\sin\alpha$ (our derivation: per unit of $x$, the edge rises $\tan\alpha$, which costs $\tan\alpha$ diagonal steps and $1 - \tan\alpha$ axis steps). That factor is 1 along the axes and the diagonals and reaches 1.082 at 22.5°, so the estimate depends on the orientation even for a perfectly digitized straight edge.

Fourier descriptors

Write the boundary as a complex function $z(t) = x(t) + i\,y(t)$ of arc length $t \in [0, L)$, sampled at $M$ points of equal arc length, $z_m = z(mL/M)$. Its Fourier descriptors are the coefficients of the discrete Fourier transform

$$c_k = \frac{1}{M} \sum_{m=0}^{M-1} z_m\, e^{-2\pi i k m / M}, \qquad z_m = \sum_{k} c_k\, e^{2\pi i k m / M}.$$

Each $c_k$ with $k \ne 0$ is a circle traversed $|k|$ times, counterclockwise for $k > 0$ and clockwise for $k < 0$; the boundary is their sum, and truncating to $|k| \le K$ gives a smoothed outline. This complex form is what OpenCV’s fourierDescriptor computes (it resamples the contour and applies a DFT), with a reference to Persoon and Fu (1977); it resamples the contour only when its length differs from the requested number of points. Belongie et al. (2002) name Zahn and Roskies (1972) and Persoon and Fu (1977) as early work that used Fourier descriptors for silhouettes.

The effect of each transformation on the coefficients follows from the DFT’s linearity and shift theorem (our derivation, which the first example checks):

Change of the contour	Effect on the coefficients
translation by $t_0$	only $c_0$ changes, by $t_0$
scaling by $s$ and rotation by $\varphi$	every $c_k$ is multiplied by $s\, e^{i\varphi}$
start point moved by $m_0$ samples, $z'_m = z_{m+m_0}$	$c_k$ is multiplied by $e^{2\pi i k m_0/M}$
reversed direction, $z'_m = z_{-m}$	$c_k$ becomes $c_{-k}$
reflection, $z \mapsto \bar z$, with the direction reversed to keep it counterclockwise	$c_k$ becomes $\bar c_k$

Dropping $c_0$ removes translation. Dividing by $|c_1|$, which for a counterclockwise contour of a compact shape is usually the largest coefficient (our reasoning), removes scale. Taking magnitudes $|c_k|/|c_1|$ removes rotation and the start point together, because both only change phases. The magnitude descriptor is cheap and needs no alignment, but it discards every phase, so it also cannot tell a shape from its mirror image (the last row of the table leaves magnitudes unchanged), and different shapes can have identical magnitudes.

Keeping the phases requires aligning them instead. OpenCV’s ContourFitting matches two descriptor sets $a_k$ and $b_k$ by searching for a scale $s$, a rotation $\varphi$ and a start-point shift $\alpha$; its documentation writes the distance as $\sum_n (a_n - s\, b_n e^{j(n\alpha + \varphi)})^2$, while its code sums the absolute values $|a_n - s\, b_n e^{j(n\alpha + \varphi)}|$ over the first harmonics and finds $\alpha$ by Newton–Raphson iterations. The examples below use a least-squares form of the same idea, our own: for each whole-sample start shift $m_0$ the best complex factor $w = s e^{i\varphi}$ has a closed form, and one FFT evaluates all shifts at once,

$$d(a, b) = \min_{m_0,\, w} \left( \frac{\sum_{1 \le |k| \le K} \left| a_k - w\, b_k\, e^{2\pi i k m_0/M} \right|^2}{\sum_{1 \le |k| \le K} |a_k|^2} \right)^{1/2}.$$

How fast the coefficients decay sets how many are needed. For a polygon parametrized by arc length, $z'$ is constant along each edge and jumps at each vertex, so $z''$ is a sum of impulses $(u_j - u_{j-1})\,\delta(t - t_j)$, where $u_j$ is the unit direction of edge $j$ and $t_j$ the arc position of vertex $j$. Differentiating the Fourier series twice gives the exact coefficients (our derivation)

$$c_k = -\frac{L}{4\pi^2 k^2} \sum_j (u_j - u_{j-1})\, e^{-2\pi i k t_j / L}, \qquad k \neq 0,$$

which fall as $1/k^2$. The mean squared error of a truncation to $|k| \le K$ is, by Parseval’s theorem, the sum of $|c_k|^2$ over the dropped terms, so the RMS error of the reconstruction falls as $K^{-3/2}$ for a shape with corners. A smooth boundary has coefficients that fall faster than any power of $k$, and few harmonics suffice.

Elliptic Fourier descriptors expand $x(t)$ and $y(t)$ separately in real Fourier series, giving four coefficients $(a_n, b_n, c_n, d_n)$ per harmonic; each harmonic traces an ellipse. They are associated with Kuhl and Giardina (1982). The pyefd library, which cites that paper, normalizes them on request (in version 1.8.0 the normalize argument of elliptic_fourier_descriptors defaults to False, although its README says the coefficients are normalized by default) by shifting the phase so that the first harmonic starts at its major axis, rotating so that this axis is parallel to the $x$ axis, forcing a counterclockwise orientation and dividing by the first harmonic’s size, after which the first three coefficients are constant and can be dropped.

Shape contexts

Belongie, Malik and Puzicha (2002) represent a shape by $n$ points sampled from its internal or external contours, for instance edge pixels from an edge detector, with roughly uniform spacing; the points need not be landmarks or curvature extrema, and they typically used about 100 (300 for the MPEG-7 silhouettes and the trademarks). The shape context of point $p_i$ is the histogram of the positions of the other $n - 1$ points relative to it (their Eq. 1),

$$h_i(k) = \#\{\, q \ne p_i : (q - p_i) \in \text{bin}(k) \,\},$$

with bins uniform in log-polar coordinates, which makes the descriptor more sensitive to nearby points than to distant ones. They used 5 bins of $\log r$ over the range 0.125 to 2 and 12 angle bins in all their experiments (their Appendix A calls the angle bins “radial”; Fig. 3 gives them as bins of $\theta$), with distances divided by the mean distance between the $n^2$ point pairs of the shape for scale invariance. Translation invariance is built in, since all positions are relative.

Two points are compared with the $\chi^2$ test statistic between their normalized histograms,

$$C_{ij} = \frac{1}{2} \sum_{k=1}^{K} \frac{\left[h_i(k) - h_j(k)\right]^2}{h_i(k) + h_j(k)},$$

the histogram distance that similarity and distance measures discusses, and two shapes are put in correspondence by the one-to-one matching $\pi$ that minimizes $H(\pi) = \sum_i C(p_i, q_{\pi(i)})$ (their Eq. 2), an assignment problem that the Hungarian algorithm solves in $O(N^3)$; the authors used a faster algorithm. They add “dummy” points with a fixed cost to absorb outliers and to make the cost matrix square when the point counts differ, and allow an appearance term in $C_{ij}$, such as the difference of tangent angles.

The full method then iterates: given correspondences, it fits a regularized thin-plate spline that maps one shape onto the other, re-matches, and after three iterations in their large-scale experiments scores the pair by a weighted sum of three terms, the symmetric shape-context matching cost, an appearance cost from grey-level windows around corresponding points, and the thin-plate spline’s bending energy (their §4–5). Recognition is nearest-neighbour search over stored prototype shapes; for COIL-20 the prototypes were chosen by k-medoids clustering of the shape distances.

Rotation invariance is optional. Measuring each point’s angles from its own tangent direction instead of the image’s $x$ axis makes the descriptor completely rotation invariant, but the authors point out that many applications need the absolute frame (distinguishing a 6 from a 9), that many points have no reliable tangent, and that appearance features lose discriminative power when measured in different frames. OpenCV’s ShapeContextDistanceExtractor, which implements the method after the paper, defaults to 12 angular and 4 radial bins between radii 0.2 and 2, three iterations, the $\chi^2$ cost and a thin-plate spline transformer.

Other families

Belongie et al. (2002) survey earlier approaches. Silhouette boundaries, a single closed curve, were described with Fourier descriptors, with the graph of the medial axis (by Kimia, Zucker and colleagues) and with dynamic-programming matching of curves by edit distance; the MPEG-7 standard’s comparison of silhouette descriptors had approaches by Latecki et al. and by Mokhtarian et al. (curvature scale space) as its leaders. Point-set methods include the Hausdorff distance of Huttenlocher et al., which can be extended to partial matching and clutter but returns no correspondences, and methods based on distance transforms; the symmetric nearest-point average used to compare point clouds is the Chamfer distance.

Properties

Invariance by construction or by normalization. The families differ in where invariance comes from:

Descriptor	Translation	Scale	Rotation	Start point	Reflection
circularity, solidity, extent, eccentricity, Euler number	yes	yes	yes, up to sampling	not applicable	yes
chain code, first difference with minimal rotation	yes	no	90° steps only	yes	no
Fourier magnitudes $\lvert c_k \rvert / \lvert c_1 \rvert$, $k \neq 0$	yes	yes	yes	yes	invariant, so mirror images are confused
Fourier descriptors with aligned distance	yes	yes	yes	yes	distinguished
shape context, absolute frame	yes	yes	no	not applicable	distinguished
shape context, tangent frame	yes	yes	yes	not applicable	distinguished

The table follows from the definitions above (our summary). The examples check the translation, scale, rotation, start-point and reflection entries of the two Fourier rows, and the translation, scale and rotation entries of the shape-context rows; the reflection entries of the shape-context rows are not tested.

Global descriptors need a clean, whole outline. A Fourier descriptor or a chain code describes one closed curve. An occluded part, a merge with a neighbouring object or a hole changes the curve and with it every coefficient, and two touching objects give one curve. Region measures have the same weakness, which is why they follow thresholding or image segmentation and a connected-component step. Shape contexts are computed and matched point by point, so, in our reading, a missing part removes some points rather than changing every descriptor completely; Belongie et al. describe them as rich descriptors that are insensitive to small perturbations of parts of the shape, state that they have no theoretical guarantees for this, and evaluated robustness to deformation, noise and outliers experimentally. (Each shape context still describes the whole shape from its point, which the authors call a globally discriminative characterization.)

Low orders describe the overall form; high orders describe corners and noise. For Fourier descriptors this is the decay derived above; it is also why a truncated series is a smoothing filter and why comparing the first few coefficients is robust to boundary noise but blind to fine detail.

Sampling is part of the descriptor. Chain codes and region perimeters are measured on the pixel grid, so they change when the shape is rotated by an angle that is not a multiple of 90°. Fourier descriptors depend on the resampling of the contour: a start point that falls between samples changes the samples themselves, not only the phases. Shape contexts depend on which points are sampled; Belongie et al. note that independent random samples on two shapes cause jitter in the correspondences, which the regularized thin-plate spline smooths.

Completeness and uniqueness. Fourier descriptors with phases are complete: with all coefficients the contour can be recovered exactly. Magnitudes are not. Region scalars are far from complete: many different shapes share a circularity and a solidity.

Examples

Chain codes and Fourier descriptors

The first script builds six closed shapes, each normalized to unit area with its centroid at the origin: a T, a U and an F (polygons; the F has no mirror symmetry), a three-lobed trefoil $r = 1 + 0.3\cos 3\theta$ and an asymmetric blob $r = 1 + 0.25\cos 2\theta + 0.15\sin 3\theta$ (smooth curves, sampled as polygons of 2,000 vertices), and a crescent (a unit disk minus a disk of radius 0.8 centred at $(0.45, 0)$, with two sharp tips). Contours run counterclockwise with $x$ to the right and $y$ up; images have row $= 96 - y$ and column $= 96 + x$, so the picture is the right way up.

For chain codes, each shape is rendered at a size of 40 pixels (the square root of its area) by testing each pixel centre, and its outer boundary is traced by Moore-neighbour tracing with Jacob’s stopping rule. The checks: the code closes; the traced pixels are exactly the object pixels with a background pixel among their four neighbours; and for rotations by one, two and three quarter turns (np.rot90) the code becomes a cyclic shift of the original plus $2q$ (mod 8) and the minimal rotation of the difference code is unchanged. The T is then rendered at sizes 20, 40 and 80 pixels and 13 angles from 0° to 90°, each with a random sub-pixel offset, and three perimeter estimates are compared with its true perimeter: the chain-code length, scikit-image’s perimeter (neighbourhood 4) and perimeter_crofton (4 directions).

For Fourier descriptors, the contours are resampled at $M$ points of equal arc length from the first vertex. The checks compare the DFT coefficients of the polygons with the closed form above, reconstruct a contour built from harmonics $-3$ to $5$ exactly with $K = 5$, and transform each shape by a scale of 2.5, a rotation of 1 radian, a translation and a start shift of 37 samples, after which the magnitude descriptor must be unchanged, the aligned distance zero and the transformation recovered. The measures: the RMS distance between each contour and its reconstruction from harmonics $|k| \le K$, for $K$ from 1 to 64; the change in the magnitude descriptor when the start point falls half a sample later, for $M$ = 64, 256 and 1024; the aligned distance between the T and the F and their mirror images, whose magnitudes are identical by the table above; and a T whose coefficients of orders 2 to 32 have random phases, so that its magnitudes equal the T’s.

The plan, written before the code ran, stated these expectations. (P1) At size 80, the chain-code error at 22.5° exceeds that at 0° by 6 to 10 percentage points, from the factor 1.082 for oblique edges and a shortening of about 4 pixels at 0°, because the chain runs through pixel centres about half a pixel inside the true boundary. (P2) At every size, the range of the Crofton estimate over the angles is smaller than that of the chain code. (FD-H1) For the T, U, F and crescent, the log-log slope of the reconstruction error between $K$ = 8 and 64 lies between −1.7 and −1.3, near the −1.5 derived above. (FD-H2) For the trefoil and the blob, it is steeper than −3, with an error below $10^{-5}$ at $K = 64$. (FD-H3) The half-sample start error falls as $M$ grows. (FD-H4) The F’s aligned distance from its mirror image is at least ten times the T’s. (FD-H5) The phase-scrambled T is at least half as far from the T as the U is.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread, set before NumPy is imported
import numpy as np
from scipy import ndimage
from skimage.measure import perimeter, perimeter_crofton

# ---------------- shapes: closed polygons, counterclockwise, x right and y up ----------------
def polar(r, n=2000):
    t = 2 * np.pi * np.arange(n) / n
    return np.c_[r(t) * np.cos(t), r(t) * np.sin(t)]

def crescent(n=1000):                         # unit disk minus a disk of radius 0.8 centred at (0.45, 0)
    a = np.arctan2(0.7806, 0.625)
    outer = np.linspace(a, 2 * np.pi - a, n)
    b = np.arctan2(0.7806, 0.175)
    inner = np.linspace(-b, -2 * np.pi + b, n)[1:-1]
    return np.r_[np.c_[np.cos(outer), np.sin(outer)], np.c_[0.45 + 0.8 * np.cos(inner), 0.8 * np.sin(inner)]]

RAW = {"T": [(1, 0), (2, 0), (2, 2), (3, 2), (3, 3), (0, 3), (0, 2), (1, 2)],
       "U": [(0, 0), (3, 0), (3, 3), (2, 3), (2, 1), (1, 1), (1, 3), (0, 3)],
       "F": [(0, 0), (1, 0), (1, 1.6), (2.2, 1.6), (2.2, 2.4), (1, 2.4), (1, 3), (3, 3), (3, 4), (0, 4)],
       "trefoil": polar(lambda t: 1 + 0.3 * np.cos(3 * t)),
       "blob": polar(lambda t: 1 + 0.25 * np.cos(2 * t) + 0.15 * np.sin(3 * t)),
       "crescent": crescent()}

def area_centroid(P):
    x, y = P[:, 0], P[:, 1]; x1, y1 = np.roll(x, -1), np.roll(y, -1)
    cr = x * y1 - x1 * y; A = cr.sum() / 2
    return A, np.array([((x + x1) * cr).sum(), ((y + y1) * cr).sum()]) / (6 * A)

def normalize(P):                             # centroid at the origin, unit area
    P = np.asarray(P, float); A, c = area_centroid(P)
    return (P - c) / np.sqrt(A)

SHAPES = {k: normalize(v) for k, v in RAW.items()}

def edges(P):
    d = np.roll(P, -1, axis=0) - P
    return d, np.hypot(d[:, 0], d[:, 1])

def resample(P, M, offset=0.0):
    """M points at equal arc length along the closed polygon, the first one offset*L/M after vertex 0."""
    d, ln = edges(P); L = ln.sum(); cum = np.r_[0, np.cumsum(ln)]
    s = ((offset + np.arange(M)) * L / M) % L
    j = np.searchsorted(cum, s, side="right") - 1
    return P[j] + d[j] * ((s - cum[j]) / ln[j])[:, None]

def simple_polygon(P):                        # no two non-adjacent edges intersect
    A, B = P, np.roll(P, -1, axis=0); n = len(P)
    def orient(p, q, r):
        return np.sign((q[..., 0] - p[..., 0]) * (r[..., 1] - p[..., 1]) - (q[..., 1] - p[..., 1]) * (r[..., 0] - p[..., 0]))
    i, j = np.triu_indices(n, 2); keep = ~((i == 0) & (j == n - 1)); i, j = i[keep], j[keep]
    hit = (orient(A[i], B[i], A[j]) != orient(A[i], B[i], B[j])) & (orient(A[j], B[j], A[i]) != orient(A[j], B[j], B[i]))
    return not hit.any()

# ---------------- Fourier descriptors ----------------
def fd(z):                                    # c_k = (1/M) sum_m z_m exp(-2 pi i k m / M)
    return np.fft.fft(z) / len(z)

def exact_fd(P, ks):
    """Exact Fourier coefficients of a polygon parametrized by arc length (period L): z'' is a sum of
    impulses (u_j - u_{j-1}) at the vertices, so c_k = -L/(4 pi^2 k^2) sum_j (u_j - u_{j-1}) exp(-2 pi i k t_j / L)."""
    d, ln = edges(P); L = ln.sum(); t = np.r_[0, np.cumsum(ln)[:-1]]
    u = (d[:, 0] + 1j * d[:, 1]) / ln; du = u - np.roll(u, 1)
    ks = np.asarray(ks, float)
    return -L / (4 * np.pi ** 2 * ks ** 2) * (du[None, :] * np.exp(-2j * np.pi * ks[:, None] * t[None, :] / L)).sum(1)

def magnitude_descriptor(c, K):              # |c_k| / |c_1| for k = -K..K, k not 0 or 1
    ks = [k for k in range(-K, K + 1) if k not in (0, 1)]
    return np.abs(c[ks]) / np.abs(c[1])

def aligned_distance(a, b, K):
    """min over start shift m0 (whole samples) and complex factor w = s e^{i phi} of sum_k |a_k - w b_k e^{2 pi i k m0/M}|^2,
    k = +-1..+-K, relative to sum |a_k|^2 (a least-squares form of the distance in OpenCV's ContourFitting)."""
    M = len(a); ks = np.r_[1:K + 1, -K:0]
    g = np.zeros(M, complex); g[ks % M] = a[ks] * np.conj(b[ks])
    cross = np.fft.fft(g)                     # cross[m0] = sum_k a_k conj(b_k) exp(-2 pi i k m0 / M)
    Ea, Eb = np.sum(np.abs(a[ks]) ** 2), np.sum(np.abs(b[ks]) ** 2)
    m0 = int(np.argmax(np.abs(cross)))
    w = cross[m0] / Eb
    res = a[ks] - w * b[ks] * np.exp(2j * np.pi * ks * m0 / M)
    return np.sqrt(np.sum(np.abs(res) ** 2) / Ea), abs(w), np.angle(w), m0

def recon_error(c, K):                        # RMS distance between the samples and their K-harmonic reconstruction
    M = len(c); k = np.fft.fftfreq(M, 1 / M)
    return np.sqrt(np.sum(np.abs(c[np.abs(k) > K]) ** 2))

# ---------------- rasterization and chain codes ----------------
N_IMG = 192
def render(P, size, angle_deg, off):          # pixel-centre test; row = N/2 - y, col = N/2 + x
    a = np.deg2rad(angle_deg); R = np.array([[np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]])
    Q = size * P @ R.T + off
    rr, cc = np.mgrid[:N_IMG, :N_IMG].astype(float)
    px, py = cc - N_IMG / 2, N_IMG / 2 - rr; inside = np.zeros((N_IMG, N_IMG), bool)
    for (x0, y0), (x1, y1) in zip(Q, np.roll(Q, -1, axis=0)):
        cross = (y0 > py) != (y1 > py)
        inside ^= cross & (px < x0 + (py - y0) * (x1 - x0) / np.where(cross, y1 - y0, 1))
    return inside

DIRS = [(0, 1), (-1, 1), (-1, 0), (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1)]   # (d_row, d_col) of codes 0..7: E, NE, N, ..., SE
CODE = {d: k for k, d in enumerate(DIRS)}
CW = [(0, -1), (-1, -1), (-1, 0), (-1, 1), (0, 1), (1, 1), (1, 0), (1, -1)]       # W, NW, N, NE, E, SE, S, SW: clockwise on screen

def trace(img):
    """Moore-neighbour tracing of the outer boundary (8-connected), with Jacob's stopping criterion.
    Returns the start pixel and the Freeman chain code."""
    r0, c0 = np.argwhere(img)[0]              # topmost, then leftmost object pixel: its west neighbour is background
    p, b = (r0, c0), (r0, c0 - 1); code, first = [], None
    while True:
        k = CW.index((b[0] - p[0], b[1] - p[1]))
        for i in range(1, 9):
            d = CW[(k + i) % 8]; q = (p[0] + d[0], p[1] + d[1])
            if img[q]:
                prev = CW[(k + i - 1) % 8]; nb = (p[0] + prev[0], p[1] + prev[1]); break
        step = (q[0] - p[0], q[1] - p[1])
        if p == (r0, c0) and first is not None and step == first:
            return (r0, c0), np.array(code)
        if first is None:
            first = step
        code.append(CODE[step]); p, b = q, nb

def decode(start, code):
    steps = np.array([DIRS[k] for k in code])
    return np.vstack([start, start + np.cumsum(steps, axis=0)])

def shape_number(code):                       # first difference, rotated to its lexicographically smallest form
    d = (code - np.roll(code, 1)) % 8
    return min(tuple(np.roll(d, -i)) for i in range(len(d)))

def is_cyclic_shift(a, b):
    return len(a) == len(b) and any(np.array_equal(np.roll(a, -i), b) for i in range(len(a)))

def chain_length(code):
    return np.sum(code % 2 == 0) + np.sqrt(2) * np.sum(code % 2 == 1)

# ======================= checks =======================
for name, P in SHAPES.items():
    A, c = area_centroid(P)
    print(f"D1 {name:8s}: {len(P):4d} vertices, simple {simple_polygon(P)}, area {A:.12f}, |centroid| {np.hypot(*c):.1e}")

rng = np.random.default_rng(0)
for name, P in SHAPES.items():
    img = render(P, 40, 0.0, rng.uniform(-0.5, 0.5, 2))
    start, code = trace(img); path = decode(start, code)
    border = img & ~ndimage.binary_erosion(img, ndimage.generate_binary_structure(2, 1), border_value=0)
    visited = np.zeros_like(img); visited[path[:, 0], path[:, 1]] = True
    n_lab = ndimage.label(img, ndimage.generate_binary_structure(2, 2))[1]
    rot = [trace(np.rot90(img, q))[1] for q in (1, 2, 3)]
    ok = [np.array_equal(path[0], path[-1]), np.array_equal(visited, border), n_lab == 1,
          all(is_cyclic_shift((code + 2 * q) % 8, r) for q, r in zip((1, 2, 3), rot)),
          all(shape_number(r) == shape_number(code) for r in rot)]
    print(f"C1 {name:8s}: {len(code)} links; closes / traced == 4-border / one component / rot90 code+2q / shape number: {ok}")

M, K = 4096, 64
ks = np.r_[1:K + 1, -K:0]
CONT = {k: (lambda z: z[:, 0] + 1j * z[:, 1])(resample(P, M)) for k, P in SHAPES.items()}
C = {k: fd(z) for k, z in CONT.items()}
print("B-F1 DFT (M=4096) vs closed-form coefficients, max |diff|/|c1| for 1<=|k|<=64 (tol 1e-5): "
      + ", ".join(f"{k} {np.max(np.abs(C[k][ks] - exact_fd(SHAPES[k], ks))) / abs(C[k][1]):.1e}" for k in ("T", "U", "F", "crescent")))
cb = np.zeros(M, complex); cb[1] = 1.0; rb = np.random.default_rng(7)
for k in (-3, -2, -1, 2, 3, 4, 5):
    cb[k % M] = 0.15 * (rb.normal() + 1j * rb.normal()) / abs(k)
cbb = fd(np.fft.ifft(cb) * M)
print(f"B-F2 band-limited contour (harmonics -3..5): error at K=5 {recon_error(cbb, 5):.1e} (tol 1e-12), at K=4 {recon_error(cbb, 4):.1e} (>= 1e-3)")
Mi, Ki = 1024, 32
Z = {k: (lambda z: z[:, 0] + 1j * z[:, 1])(resample(P, Mi)) for k, P in SHAPES.items()}
Ci = {k: fd(z) for k, z in Z.items()}
worst = [0.0, 0.0, 0.0, True]
for name, z in Z.items():
    a = fd(2.5 * np.exp(1j * 1.0) * np.roll(z, -37) + (3 - 2j))
    d, s, ph, m0 = aligned_distance(a, Ci[name], Ki)
    worst = [max(worst[0], np.max(np.abs(magnitude_descriptor(a, Ki) - magnitude_descriptor(Ci[name], Ki)))),
             max(worst[1], d), max(worst[2], abs(s - 2.5), abs(ph - 1.0)), worst[3] and m0 == 37]
print(f"B-F3 copy with s=2.5, phi=1 rad, t=3-2i, start +37 samples (all shapes): magnitude max diff {worst[0]:.1e} (tol 1e-12), "
      f"aligned distance {worst[1]:.1e} (tol 1e-12), s/phi error {worst[2]:.1e} (tol 1e-10), m0 = 37 {worst[3]}")
mirror = {k: fd(np.conj(Z[k][::-1])) for k in ("T", "F")}
print("B-F4 mirror image, magnitude max diff (tol 1e-12): "
      + ", ".join(f"{k} {np.max(np.abs(magnitude_descriptor(mirror[k], Ki) - magnitude_descriptor(Ci[k], Ki))):.1e}" for k in mirror))
scr = Ci["T"].copy(); psi = np.random.default_rng(5).uniform(0, 2 * np.pi, 2 * Ki + 1); kk = np.fft.fftfreq(Mi, 1 / Mi)
sel = (np.abs(kk) >= 2) & (np.abs(kk) <= Ki); scr[sel] *= np.exp(1j * psi[(kk[sel] + Ki).astype(int)])
print(f"B-F5 phase-scrambled T, magnitude max diff from T {np.max(np.abs(magnitude_descriptor(scr, Ki) - magnitude_descriptor(Ci['T'], Ki))):.1e} (tol 1e-12)")

# ======================= chain-code perimeter =======================
print("\nT shape: perimeter estimates, relative error (%) against the true perimeter; chain = n_even + sqrt(2) n_odd,")
print("skimage = measure.perimeter(img, 4), Crofton = measure.perimeter_crofton(img, 4)")
ANG = np.arange(13) * 7.5; SIZES = [20, 40, 80]; L_T = edges(SHAPES["T"])[1].sum()
offs = np.random.default_rng(1).uniform(-0.5, 0.5, (len(SIZES), len(ANG), 2)); res = {}
for i, s in enumerate(SIZES):
    for j, a in enumerate(ANG):
        img = render(SHAPES["T"], s, a, offs[i, j]); code = trace(img)[1]; true = L_T * s
        res[s, a] = [100 * (chain_length(code) / true - 1), 100 * (perimeter(img, 4) / true - 1),
                     100 * (perimeter_crofton(img, 4) / true - 1)]
print("angle " + "".join(f"  chain{s:<3d} skimg{s:<3d} crof{s:<3d}" for s in SIZES))
for a in ANG:
    print(f"{a:5.1f} " + "".join(f"{res[s, a][0]:9.2f}{res[s, a][1]:9.2f}{res[s, a][2]:8.2f}" for s in SIZES))
for s in SIZES:
    rg = [max(res[s, a][m] for a in ANG) - min(res[s, a][m] for a in ANG) for m in range(3)]
    print(f"size {s}: range over angles (points): chain {rg[0]:.2f}, skimage {rg[1]:.2f}, Crofton {rg[2]:.2f}; "
          f"chain at 22.5 minus at 0: {res[s, 22.5][0] - res[s, 0.0][0]:.2f}")

# ======================= Fourier descriptors =======================
KS = [1, 2, 4, 8, 16, 32, 64]
print("\nReconstruction from harmonics |k| <= K: RMS point error in units of sqrt(area) (M = 4096)")
print("shape    " + "".join(f"{'K=' + str(k):>10s}" for k in KS) + "   slope K=8..64")
for name in SHAPES:
    err = [recon_error(C[name], k) for k in KS]
    print(f"{name:8s} " + "".join(f"{e:10.2e}" for e in err) + f"   {np.polyfit(np.log(KS[3:]), np.log(err[3:]), 1)[0]:.3f}")

print("\nStart point between samples: max |difference| of the magnitude descriptor (K = 16), offset 0.5 vs 0 samples")
for name in ("T", "trefoil"):
    out = []
    for Mm in (64, 256, 1024):
        z0, z5 = [(lambda z: z[:, 0] + 1j * z[:, 1])(resample(SHAPES[name], Mm, o)) for o in (0.0, 0.5)]
        out.append(np.max(np.abs(magnitude_descriptor(fd(z5), 16) - magnitude_descriptor(fd(z0), 16))))
    print(f"  {name:8s} M=64 {out[0]:.2e}   M=256 {out[1]:.2e}   M=1024 {out[2]:.2e}")

print("\nAligned distance (start, rotation, scale removed; M = 1024, K = 32)")
print("  mirror image: " + ", ".join(f"{k} {aligned_distance(mirror[k], Ci[k], Ki)[0]:.2e}" for k in mirror))
print("  from T: " + ", ".join(f"{k} {aligned_distance(Ci[k], Ci['T'], Ki)[0]:.3f}" for k in SHAPES if k != "T")
      + f", phase-scrambled T {aligned_distance(scr, Ci['T'], Ki)[0]:.3f}")

D1 T       :    8 vertices, simple True, area 1.000000000000, |centroid| 5.6e-17
D1 U       :    8 vertices, simple True, area 1.000000000000, |centroid| 4.6e-17
D1 F       :   10 vertices, simple True, area 1.000000000000, |centroid| 8.9e-17
D1 trefoil : 2000 vertices, simple True, area 1.000000000000, |centroid| 8.3e-17
D1 blob    : 2000 vertices, simple True, area 1.000000000000, |centroid| 1.9e-16
D1 crescent: 1998 vertices, simple True, area 1.000000000000, |centroid| 1.3e-16
C1 T       : 208 links; closes / traced == 4-border / one component / rot90 code+2q / shape number: [True, True, True, True, True]
C1 U       : 238 links; closes / traced == 4-border / one component / rot90 code+2q / shape number: [True, True, True, True, True]
C1 F       : 243 links; closes / traced == 4-border / one component / rot90 code+2q / shape number: [True, True, True, True, True]
C1 trefoil : 145 links; closes / traced == 4-border / one component / rot90 code+2q / shape number: [True, True, True, True, True]
C1 blob    : 136 links; closes / traced == 4-border / one component / rot90 code+2q / shape number: [True, True, True, True, True]
C1 crescent: 216 links; closes / traced == 4-border / one component / rot90 code+2q / shape number: [True, True, True, True, True]
B-F1 DFT (M=4096) vs closed-form coefficients, max |diff|/|c1| for 1<=|k|<=64 (tol 1e-5): T 2.0e-07, U 4.8e-07, F 3.5e-07, crescent 1.5e-07
B-F2 band-limited contour (harmonics -3..5): error at K=5 1.8e-16 (tol 1e-12), at K=4 2.8e-02 (>= 1e-3)
B-F3 copy with s=2.5, phi=1 rad, t=3-2i, start +37 samples (all shapes): magnitude max diff 2.8e-16 (tol 1e-12), aligned distance 3.6e-16 (tol 1e-12), s/phi error 8.9e-16 (tol 1e-10), m0 = 37 True
B-F4 mirror image, magnitude max diff (tol 1e-12): T 2.8e-17, F 1.7e-16
B-F5 phase-scrambled T, magnitude max diff from T 5.6e-17 (tol 1e-12)

T shape: perimeter estimates, relative error (%) against the true perimeter; chain = n_even + sqrt(2) n_odd,
skimage = measure.perimeter(img, 4), Crofton = measure.perimeter_crofton(img, 4)
angle   chain20  skimg20  crof20   chain40  skimg40  crof40   chain80  skimg80  crof80 
  0.0     -4.20    -4.20   -6.67    -2.72    -2.72   -6.52    -1.51    -1.51   -6.00
  7.5     -2.75    -2.75   -5.30     2.03     2.03   -2.02     3.09     3.05   -1.64
 15.0      1.34     1.34   -1.43     3.98     3.98   -0.18     6.30     6.11    1.41
 22.5      3.56     2.21    0.68     6.20     5.81    1.93     6.94     6.94    2.01
 30.0      3.92     2.76    1.02     4.85     4.66    0.66     5.70     5.70    0.84
 37.5     -0.45    -0.83   -3.12     1.92     1.72   -2.13     3.43     3.33   -1.32
 45.0     -4.04    -4.04   -6.53    -1.95    -1.95   -5.79    -0.91    -0.91   -5.43
 52.5      0.87     0.49   -1.87     1.53     1.34   -2.49     3.76     3.66   -1.01
 60.0      3.92     3.15    1.02     4.08     4.08   -0.08     6.17     6.07    1.28
 67.5      0.93     0.93   -1.82     5.92     5.35    1.67     6.80     6.80    1.88
 75.0      0.25     0.25   -2.46     6.22     5.55    1.96     5.37     5.37    0.52
 82.5      1.75     1.56   -1.03     2.03     1.93   -2.02     3.56     3.51   -1.19
 90.0     -4.20    -4.20   -6.67    -1.79    -1.79   -5.64    -0.58    -0.58   -5.12
size 20: range over angles (points): chain 8.12, skimage 7.34, Crofton 7.69; chain at 22.5 minus at 0: 7.76
size 40: range over angles (points): chain 8.94, skimage 8.53, Crofton 8.48; chain at 22.5 minus at 0: 8.91
size 80: range over angles (points): chain 8.45, skimage 8.45, Crofton 8.01; chain at 22.5 minus at 0: 8.45

Reconstruction from harmonics |k| <= K: RMS point error in units of sqrt(area) (M = 4096)
shape           K=1       K=2       K=4       K=8      K=16      K=32      K=64   slope K=8..64
T          2.06e-01  1.22e-01  4.79e-02  1.60e-02  6.69e-03  2.28e-03  8.58e-04   -1.423
U          3.37e-01  1.40e-01  6.10e-02  2.04e-02  7.63e-03  2.75e-03  9.78e-04   -1.463
F          2.77e-01  1.59e-01  7.43e-02  2.06e-02  9.11e-03  3.24e-03  1.10e-03   -1.417
trefoil    1.15e-01  5.99e-02  1.45e-02  3.26e-03  6.83e-04  3.52e-05  7.86e-07   -4.033
blob       7.01e-02  4.92e-02  6.72e-03  2.25e-03  3.78e-04  2.98e-05  6.45e-07   -3.896
crescent   3.12e-01  1.08e-01  4.19e-02  1.34e-02  5.14e-03  1.90e-03  6.69e-04   -1.441

Start point between samples: max |difference| of the magnitude descriptor (K = 16), offset 0.5 vs 0 samples
  T        M=64 7.49e-04   M=256 3.87e-05   M=1024 2.32e-06
  trefoil  M=64 3.52e-07   M=256 9.30e-08   M=1024 1.13e-08

Aligned distance (start, rotation, scale removed; M = 1024, K = 32)
  mirror image: T 2.22e-03, F 3.96e-01
  from T: U 0.547, F 0.565, trefoil 0.217, blob 0.285, crescent 0.479, phase-scrambled T 0.299

Every check passes: the tracer, the difference-code normalization, the closed-form coefficients (agreement within $5 \times 10^{-7}$ of $|c_1|$; in our reading, the remaining difference is aliasing from the finite $M$) and the exact invariance of the Fourier descriptors to translation, scale, rotation and whole-sample start shifts, to rounding error.

Perimeter (P1 holds; P2 holds, narrowly). No estimator is independent of orientation. The chain-code error rises from −1.51% at 0° to 6.94% at 22.5° at size 80, a difference of 8.45 points, within the planned 6 to 10, and the pattern repeats around 67.5°; at 0°, 45° and 90° all three estimators are at their lowest. Crofton’s range over angles is smaller than the chain code’s at every size, but by less than half a point each time (7.69 against 8.12 points at size 20, 8.48 against 8.94 at 40, 8.01 against 8.45 at 80), so none of the three removes the orientation dependence. Crofton’s estimate is also biased low: −6.67%, −6.52% and −6.00% at 0° for the three sizes, and below the truth at 26 of the 39 combinations of angle and size (our count); the run does not show the source of the bias. scikit-image’s perimeter is close to the chain code at most angles and identical at 0°, 45° and 90°. A circularity computed from any of them inherits twice the relative error, so a T rotated by 22.5° appears about 15% less circular than the same T at 0° when the chain-code length is used (our arithmetic from the size-80 rows: $(0.9849/1.0694)^2 = 0.85$).

Reconstruction (FD-H1 and FD-H2 hold). The four shapes with corners have slopes between −1.463 and −1.417, close to the −1.5 of the $1/k^2$ decay; the error is still 0.09% of the shape’s size at $K = 64$ for the T. The two smooth shapes have slopes of −4.033 and −3.896 and errors below $10^{-6}$ at $K = 64$; from $K = 8$ on their error falls ever faster as $K$ grows, by a factor of 13 to 19 between $K = 16$ and 32 and of about 45 between 32 and 64. Among these six shapes, corners, not overall complexity, decided how many descriptors a shape needs.

Start point between samples (FD-H3 holds). A start point half a sample later changes the T’s magnitude descriptor by $7.5 \times 10^{-4}$ with 64 samples and $2.3 \times 10^{-6}$ with 1,024; for the trefoil the changes are 200 to 2,100 times smaller. In our reading, this follows from the resampling: shifting the samples changes which points sit near the corners, and only a smooth curve is sampled almost equally well from every start.

Mirror images and phases (FD-H4 holds; FD-H5 holds, narrowly). The F and its mirror image have identical magnitude descriptors (check B-F4) but an aligned distance of 0.396, while the T, which has a mirror axis, is 0.0022 from its mirror, only the resampling difference. The phase-scrambled T has exactly the T’s magnitudes (check B-F5) and an aligned distance of 0.299 from it, above the planned half of the T–U distance (0.274) but only just; it is farther from the T than the trefoil (0.217) and the blob (0.285) are. A magnitude-only descriptor would call it a perfect match.

The code uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2, with one BLAS thread. It ran in about 1.5 seconds; two runs printed identical output.

Shape contexts under deformation and rotation

The second script implements a reduced version of Belongie et al.’s method: the shape context with their bins (5 log-radius bins from 0.125 to 2 times the mean distance between point pairs, here the mean over the $n(n-1)$ pairs of distinct points rather than the paper’s $n^2$ pairs, and 12 angle bins; points outside the radial range are not counted), the $\chi^2$ cost and the optimal one-to-one assignment by SciPy’s linear_sum_assignment. It omits the thin-plate spline iterations, the dummy points and the appearance term, and scores a pair by the mean $\chi^2$ cost over the assignment, which ranks prototypes in the same order as the summed cost that the authors used, after one iteration with no transformation, in their rotation experiment (their Appendix A).

Each of the six shapes from the first example has a prototype of $n = 100$ points at equal arc length. Queries are drawn as 100 points at equal arc length with a random start, moved by a smooth random warp (three sinusoids per coordinate, 0.3 to 0.8 cycles per unit length, random directions and phases) whose RMS displacement is $A$ times the square root of the area, plus jitter of 0.005, then scaled by a random factor between 0.5 and 2 and translated. Because each query point is a displaced copy of a known point on the undeformed shape, the correct correspondence is known: the correspondence error is the distance between the prototype point a query point is matched to and the query point’s undeformed source, in units of the square root of the area, the median over points and then over queries. The script also prints it in units of the shape’s sample spacing.

Each query is classified by the prototype with the lowest matching cost. Ten queries per shape are drawn at $A$ = 0, 0.05, 0.1 and 0.2. A second set of ten per shape at $A = 0.05$ is rotated about its centroid by 0°, 45°, 90° and 180° and classified with absolute angles and with angles measured from each point’s tangent, estimated from its two neighbours along the contour. The checks come first. They run on prototypes with a small jitter, because exact samples of axis-parallel edges put point pairs exactly on the boundary between angle bins (at 0° and 90°), where floating-point rounding decides the bin.

The plan’s expectations: (SC-H1) the median cost of a query to its own prototype rises at each step of $A$; (SC-H2) all 60 queries are classified correctly at $A$ = 0.05 and 0.1, and accuracy at 0.2 is not above that at 0.1; (SC-H3) the median correspondence error rises at each step of $A$ and stays below one sample spacing (at least 0.0383) at $A = 0.05$; (SC-H4) with absolute angles, accuracy at 90° and at 180° is below that at 0°; (SC-H5) with tangent angles, every cost is identical at all four angles; (SC-H6) at 0°, the median margin between the nearest other prototype and the own prototype is smaller with tangent angles than with absolute ones.

The script starts with the definitions of the first one, everything above its # ======================= checks ======================= line, and continues:

# ---------------- shape contexts (Belongie, Malik and Puzicha, 2002) ----------------
from itertools import permutations
from scipy.optimize import linear_sum_assignment
NR, NT, R_IN, R_OUT = 5, 12, 0.125, 2.0      # 5 log-r bins over [0.125, 2] x mean distance, 12 angle bins

def shape_context(Q, relative=False):
    """Log-polar histograms, one per point, of the positions of the other points, normalized to sum 1.
    relative=True measures angles from each point's tangent (from its two contour neighbours)."""
    n = len(Q); D = Q[None, :, :] - Q[:, None, :]           # D[i, j] = q_j - q_i
    r = np.hypot(D[..., 0], D[..., 1]); off = ~np.eye(n, dtype=bool)
    r = r / r[off].mean()                                  # scale: mean distance over the n(n-1) pairs
    th = np.arctan2(D[..., 1], D[..., 0])
    if relative:
        tg = np.roll(Q, -1, axis=0) - np.roll(Q, 1, axis=0)
        th = th - np.arctan2(tg[:, 1], tg[:, 0])[:, None]
    rb = np.floor((np.log(np.where(off, r, 1)) - np.log(R_IN)) / (np.log(R_OUT / R_IN) / NR)).astype(int)
    tb = np.floor((th % (2 * np.pi)) / (2 * np.pi / NT)).astype(int) % NT
    ok = off & (rb >= 0) & (rb < NR)
    H = np.zeros((n, NR * NT))
    rows = np.repeat(np.arange(n)[:, None], n, 1)
    np.add.at(H, (rows[ok], (rb * NT + tb)[ok]), 1)
    return H / H.sum(1, keepdims=True)

def chi2_costs(H, G):                         # C_ij = 1/2 sum_k (h_ik - g_jk)^2 / (h_ik + g_jk), 0/0 = 0
    num = (H[:, None, :] - G[None, :, :]) ** 2; den = H[:, None, :] + G[None, :, :]
    return 0.5 * np.sum(np.where(den > 0, num / np.where(den > 0, den, 1), 0), axis=2)

def match(HQ, HP):                            # optimal one-to-one assignment; mean cost per point
    C = chi2_costs(HQ, HP); i, j = linear_sum_assignment(C)
    return C[i, j].mean(), j

def warp_field(rng, A):
    """Smooth random displacement: per coordinate, 3 sinusoids of random direction, 0.3-0.8 cycles per unit
    length and random phase, scaled so that the RMS displacement is A (units of sqrt(area))."""
    f = rng.uniform(0.3, 0.8, (2, 3))[..., None] * np.stack([np.cos(t := rng.uniform(0, 2 * np.pi, (2, 3))), np.sin(t)], -1)
    ph = rng.uniform(0, 2 * np.pi, (2, 3))
    def D(S):
        arg = 2 * np.pi * np.einsum("cjd,nd->cnj", f, S) + ph[:, None, :]
        return A * np.sqrt(2.0 / 3.0) * np.sin(arg).sum(-1).T / np.sqrt(2)
    return D

def make_query(P, n, rng, A, jitter=0.005):
    """n points at equal arc length with a random start, displaced by a smooth warp of RMS size A plus jitter,
    then scaled and translated at random. Returns the points and their undeformed sources on P."""
    S = resample(P, n, rng.uniform(0, 1))
    Q = S + warp_field(rng, A)(S) + rng.normal(0, jitter, S.shape)
    return rng.uniform(0.5, 2.0) * Q + rng.uniform(-5, 5, 2), S

# ======================= checks =======================
n = 100; names = list(SHAPES)
PROTO = {k: resample(P, n) for k, P in SHAPES.items()}
SPACING = {k: edges(P)[1].sum() / n for k, P in SHAPES.items()}
HP = {k: shape_context(Q) for k, Q in PROTO.items()}
HPr = {k: shape_context(Q, True) for k, Q in PROTO.items()}
rng = np.random.default_rng(3); b = [0.0, 0.0, True, 0.0, 1.0]
a73 = np.deg2rad(73.0); R73 = np.array([[np.cos(a73), -np.sin(a73)], [np.sin(a73), np.cos(a73)]])
for name, Q0 in PROTO.items():                  # generic copies: exact samples of straight edges sit on bin edges
    Q = Q0 + rng.normal(0, 0.005, Q0.shape); H0, H0r = shape_context(Q), shape_context(Q, True)
    Q2 = 1.7 * Q + np.array([4.0, -2.5]); perm = rng.permutation(n)
    cost, j = match(shape_context(Q2[perm]), H0)
    b = [max(b[0], np.abs(shape_context(Q2) - H0).max()), max(b[1], cost), b[2] and np.array_equal(j, perm),
         max(b[3], np.abs(shape_context(Q @ R73.T, True) - H0r).max()), min(b[4], np.abs(shape_context(Q @ R73.T) - H0).max())]
print(f"B-S1 scale and translation: max |dH| {b[0]:.1e} (tol 1e-12)")
print(f"B-S2 permuted copy: matching cost {b[1]:.1e} (tol 1e-12), permutation recovered {b[2]}")
print(f"B-S3 rotation by 73 deg: relative frame max |dH| {b[3]:.1e} (tol 1e-12); absolute frame min over shapes of max |dH| {b[4]:.2f} (>= 0.1)")
C7 = np.random.default_rng(4).random((7, 7)); i7, j7 = linear_sum_assignment(C7)
bf = min(sum(C7[i, p[i]] for i in range(7)) for p in permutations(range(7)))
print(f"B-S4 linear_sum_assignment vs brute force, 7x7: {abs(C7[i7, j7].sum() - bf):.1e} (tol 1e-12)")
rw = np.random.default_rng(12); Sx = resample(SHAPES["T"], n)
print("D-S6 warp RMS / A, mean over 200 fields (in [0.9, 1.1]): " + ", ".join(
    f"A={A}: {np.mean([np.sqrt(np.mean(np.sum(warp_field(rw, A)(Sx) ** 2, 1))) / A for _ in range(200)]):.3f}" for A in (0.05, 0.1, 0.2)))

def classify(Q, S, name, H):
    HQ = shape_context(Q, H is HPr)
    res = {k: match(HQ, H[k]) for k in names}
    own = res[name][0]; other = min(res[k][0] for k in names if k != name)
    corr = np.median(np.hypot(*(PROTO[name][res[name][1]] - S).T))
    return own < other, own, other, corr, np.array([res[k][0] for k in names])

# ======================= deformation =======================
rng = np.random.default_rng(11); AS = [0.0, 0.05, 0.1, 0.2]; out = {}
for A in AS:
    out[A] = [classify(*make_query(SHAPES[name], n, rng, A), name, HP) + (SPACING[name],) for name in names for q in range(10)]
r0 = out[0.0]
print(f"B-S5 A=0: accuracy {sum(r[0] for r in r0)}/60 (60/60), median correspondence error {np.median([r[3] for r in r0]):.4f} (<= 0.019)")
print("\nDeformation (absolute frame, 10 queries per shape, 60 in all)")
print("    A  correct  own cost  nearest other  corr. error  corr. error / spacing   (medians over queries)")
for A in AS:
    r = out[A]
    print(f"{A:5.2f}  {sum(x[0] for x in r):3d}/60  {np.median([x[1] for x in r]):8.4f}  {np.median([x[2] for x in r]):13.4f}"
          f"  {np.median([x[3] for x in r]):11.4f}  {np.median([x[3] / x[5] for x in r]):21.2f}")
print("errors at A=0.2 (true -> chosen): " + (", ".join(
    f"{names[i // 10]}->{names[int(np.argmin(out[0.2][i][4]))]}" for i in range(60) if not out[0.2][i][0]) or "none"))

# ======================= rotation =======================
rng = np.random.default_rng(13); base = [(name, *make_query(SHAPES[name], n, rng, 0.05)) for name in names for q in range(10)]
print("\nRotation (A = 0.05, the same 60 queries rotated about their centroid)")
print("frame     angle  correct  own cost  margin (nearest other - own)")
costs = {}
for label, H in (("absolute", HP), ("relative", HPr)):
    for ang in (0, 45, 90, 180):
        t = np.deg2rad(ang); R = np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])
        r = []
        for name, Q, S in base:
            c = Q.mean(0); r.append(classify((Q - c) @ R.T + c, S, name, H))
        costs[label, ang] = np.array([x[4] for x in r])
        print(f"{label:8s} {ang:6d}  {sum(x[0] for x in r):3d}/60  {np.median([x[1] for x in r]):8.4f}  {np.median([x[2] - x[1] for x in r]):8.4f}")
print(f"B-S7 relative frame, max |cost(angle) - cost(0)| over queries and prototypes: "
      f"{max(np.abs(costs['relative', a] - costs['relative', 0]).max() for a in (45, 90, 180)):.1e} (tol 1e-12)")

B-S1 scale and translation: max |dH| 0.0e+00 (tol 1e-12)
B-S2 permuted copy: matching cost 0.0e+00 (tol 1e-12), permutation recovered True
B-S3 rotation by 73 deg: relative frame max |dH| 0.0e+00 (tol 1e-12); absolute frame min over shapes of max |dH| 0.24 (>= 0.1)
B-S4 linear_sum_assignment vs brute force, 7x7: 0.0e+00 (tol 1e-12)
D-S6 warp RMS / A, mean over 200 fields (in [0.9, 1.1]): A=0.05: 0.994, A=0.1: 0.977, A=0.2: 0.978
B-S5 A=0: accuracy 60/60 (60/60), median correspondence error 0.0113 (<= 0.019)

Deformation (absolute frame, 10 queries per shape, 60 in all)
    A  correct  own cost  nearest other  corr. error  corr. error / spacing   (medians over queries)
 0.00   60/60    0.0346         0.2816       0.0113                   0.22
 0.05   60/60    0.0686         0.2815       0.0134                   0.27
 0.10   60/60    0.1389         0.2936       0.0195                   0.41
 0.20   56/60    0.2417         0.3245       0.0685                   1.33
errors at A=0.2 (true -> chosen): trefoil->blob, trefoil->F, trefoil->blob, blob->trefoil

Rotation (A = 0.05, the same 60 queries rotated about their centroid)
frame     angle  correct  own cost  margin (nearest other - own)
absolute      0   60/60    0.0724    0.1999
absolute     45   18/60    0.3610   -0.0608
absolute     90   20/60    0.3497   -0.1335
absolute    180   18/60    0.3650   -0.0797
relative      0   60/60    0.0872    0.1170
relative     45   60/60    0.0872    0.1170
relative     90   60/60    0.0872    0.1170
relative    180   60/60    0.0872    0.1170
B-S7 relative frame, max |cost(angle) - cost(0)| over queries and prototypes: 0.0e+00 (tol 1e-12)

The checks pass: the descriptor is unchanged by scale and translation, a permuted copy is matched back to its exact permutation at zero cost, the tangent frame is unchanged by rotation and the absolute frame is not, SciPy’s assignment matches brute force, and the warps have the planned size.

Deformation (SC-H1, SC-H2 and SC-H3 hold). The median cost of a query to its own prototype about doubles at each of the first two steps of $A$ and rises by three quarters at the last, from 0.035 to 0.242, while the median cost to the nearest wrong prototype moves only from 0.282 to 0.325, so the gap between the two medians shrinks from about 0.25 to about 0.08. Accuracy is 60 of 60 up to $A = 0.1$ and 56 of 60 at $A = 0.2$; of the four errors, two trefoils were taken for the blob, one trefoil for the F and one blob for the trefoil, so three of the four confuse the two smooth, rounded shapes. Correspondences are good while the shapes are recognizable: at $A = 0.05$ a query point is matched a median of 0.27 sample spacings from its true partner, and even with no warp the error is 0.22 spacings, which in our reading is the offset between two independent samplings of the same curve. At $A = 0.2$ it rises to 1.33 spacings. Without the thin-plate spline step, nothing pulls the correspondences towards a smooth mapping; in the full method, in our reading, the warp estimated from the first matches would reduce both the cost and the correspondence error of the deformed queries.

Rotation (SC-H4, SC-H5 and SC-H6 hold). With absolute angles, accuracy falls from 60 of 60 at 0° to 18 of 60 at 45°, 20 of 60 at 90° and 18 of 60 at 180°, and the median margin turns negative: a rotated shape looks more like some other upright shape than like itself. With tangent angles, every cost is the same at all four angles to the last bit (check B-S7), and all 60 queries are classified correctly at every angle. That invariance has a price at 0°: the median margin is 0.117 with tangent angles against 0.200 with absolute ones, and the own-prototype cost is higher (0.087 against 0.072), so the shapes are less well separated when nothing needs rotation invariance, which is the trade-off Belongie et al. describe.

The code uses NumPy 2.3.5 and SciPy 1.16.3, with one BLAS thread. It ran in about 14 seconds; two runs printed identical output. Both examples ran as piloted, with no change after the first run.

Common Misconceptions

• “Circularity is a precise number for a digital region.” It depends on a perimeter estimate, and every estimator in the first example changed by 7 to 9 percentage points with the T’s orientation alone. Comparing circularities of objects at different orientations, or estimated with different functions, compares the estimators as much as the shapes.
• “Fourier descriptor magnitudes identify a shape.” They are invariant to rotation and start point because they discard all phases, and with the phases goes the ability to tell a shape from its mirror image or from other shapes with the same spectrum; the phase-scrambled T in the example has exactly the T’s magnitudes. Keeping the phases and aligning them, as OpenCV’s ContourFitting does, restores the distinction.
• “More harmonics are always needed for a faithful outline.” Smooth outlines are reproduced to within $10^{-6}$ of their size by 64 harmonics on each side; the outlines that need many are those with corners, whose coefficients fall only as $1/k^2$.
• “Shape contexts are rotation invariant.” In their standard form they use absolute angles and are not; in the example, rotating the queries by 45°, 90° or 180° cut accuracy to a third or less. Measuring angles from the tangent makes them invariant but also less discriminative, and Belongie et al. recommend it only where the application needs it.
• “A chain code is a rotation-invariant description because its first difference is.” The first difference is invariant to rotations of the pixel grid by multiples of 90°. A shape rotated by any other angle is digitized into a different set of pixels, with a different code.
• “Shape descriptors are moment invariants.” Moments are one family of region descriptors, with their own invariance theory (Image Moments); contour and point-set descriptors work on the boundary or on edge points, and the shape context describes each point separately.

Where It Is Used

• Measuring segmented objects. In industrial visual inspection, blob analysis accepts or rejects each thresholded blob by measures such as its area and length; scikit-image’s regionprops computes the region descriptors above for every labelled object.
• Features for classifiers. The semi-supervised learning example classifies shapes from solidity, circularity, eccentricity and logarithms of Hu moments, and the two-stage model in the open-vocabulary segmentation example adds extent and compactness to the features of each proposal.
• Character and digit recognition. Shape contexts with a thin-plate spline gave Belongie et al. a 0.63% test error on MNIST, with 20,000 training digits and a 3-nearest-neighbour classifier, below the 0.7% they cite as the lowest then published (a boosted LeNet-4 trained with synthetic distortions); a footnote adds that DeCoste and Schölkopf reported 0.56% with virtual support vectors and the full training set.
• Silhouette and trademark retrieval. On the MPEG-7 silhouette database (Core Experiment CE-Shape-1, part B: 1,400 images of 70 classes, scored by the “bullseye” test, which counts, for each image used as a query, the correct images among its 40 best matches) Belongie et al., with 300 sample points and the lowest distance over a reference shape and its two flipped copies, reported 76.51%, against 76.45% for Latecki et al. and 75.44% for Mokhtarian et al., the best published results at the time. They also retrieved trademarks from a set of 300 with an affine transformation model.
• 3D objects from views. On COIL-20 they reported a 2.4% error with an average of only four prototype views per object chosen by k-medoids (§6.2); their introduction gives the same result as 2.5%.
• Descriptors that borrowed the log-polar layout. Mikolajczyk and Schmid’s comparison of local feature descriptors included shape context, and Dalal and Triggs tested log-polar circular blocks, reminiscent of shape contexts, in histograms of oriented gradients. In pictorial structures, Andriluka, Roth and Schiele built each part detector on densely sampled shape context descriptors.