Image Moments

Image moments are weighted sums of powers of the pixel coordinates over an image or a region. The lowest orders are familiar quantities: the zeroth moment of a binary region is its area, the first moments give its centroid, and the second moments give the orientation and elongation of the ellipse that best matches it. Combinations of higher moments can be made unchanged by translation, rotation and scaling, which made moment invariants an early shape feature for recognizing objects regardless of their pose. This article defines raw, central and normalized moments, Hu’s seven invariants, complex moments and Zernike moments, and measures how far the invariants actually hold on pixels, where sampling, thresholding and noise break the exact invariance that holds for continuous shapes.

The name comes from mechanics and statistics. Read as a mass distribution, an image’s moments are the mass, the centre of mass and the moments of inertia; read as a probability distribution over pixel positions (after dividing by the total), they are the mean and the covariance of position. Both readings are used in practice.

Definition

For an image $f(x, y)$, the raw (or geometric) moment of order $p + q$ is the sum of $x^p y^q$ weighted by the image values. The image can be a binary mask of a region, in which case the moments describe its shape, or a grey-level image or probability map, in which case they also depend on intensity. Central moments are the same sums taken about the centroid, so they do not change when the region moves; normalized central moments are divided by a power of the area, so they do not change when it is scaled. Moment invariants are functions of moments that are also unchanged by rotation, such as the seven invariants of Hu (1962). Orthogonal moments, such as Zernike moments, project the image onto a family of orthogonal polynomials instead of plain powers, which makes each moment carry independent information and lets the image be reconstructed from them.

Intuition

Imagine cutting a region out of cardboard of uniform thickness. Its weight is the zeroth moment. The point where it balances on a pin is the centroid, the first moments divided by the zeroth. Spinning it about an axis through the centroid takes more effort the farther the material lies from the axis; that resistance, the moment of inertia, is a second moment, and it is smallest about the region’s long axis. Second moments therefore say which way a region points and how elongated it is, which is why tracking, the deslanting of handwritten digits and oriented boxes use them.

Higher moments describe finer structure. Third moments measure lopsidedness (skewness): a region with more material on one side of the centroid than the other. Fourth and higher moments respond to ever smaller details of the boundary. None of these numbers alone is unchanged when the region rotates, since the coordinates themselves rotate; invariants are built by combining moments so that the rotation cancels, much as the length of a vector is unchanged when the vector turns.

Formal Definition

Raw, central and normalized moments

For an image $f$ on a pixel grid, with $x$ the column and $y$ the row,

$$m_{pq} = \sum_{x}\sum_{y} x^p\, y^q\, f(x, y), \qquad p, q = 0, 1, 2, \dots$$

and for a continuous image the sums become integrals. The area (or total mass) is $m_{00}$ and the centroid is $(\bar x, \bar y) = (m_{10}/m_{00},\, m_{01}/m_{00})$. Central moments are

$$\mu_{pq} = \sum_{x}\sum_{y} (x - \bar x)^p\, (y - \bar y)^q\, f(x, y),$$

so $\mu_{10} = \mu_{01} = 0$ and $\mu_{00} = m_{00}$. Scaling the coordinates by $s$ multiplies $\mu_{pq}$ by $s^{p+q+2}$ and $\mu_{00}$ by $s^2$, so the normalized central moments

$$\eta_{pq} = \frac{\mu_{pq}}{\mu_{00}^{\,1 + (p+q)/2}}, \qquad p + q \ge 2,$$

are unchanged by translation and uniform scaling. This is the normalization in scikit-image’s moments_normalized and in OpenCV’s cv::moments (fields nu20 to nu03).

Second moments: orientation and the equivalent ellipse

Divided by $\mu_{00}$, the second central moments form the covariance matrix of pixel position,

$$\mathbf C = \frac{1}{\mu_{00}} \begin{pmatrix} \mu_{20} & \mu_{11} \\ \mu_{11} & \mu_{02} \end{pmatrix}.$$

Its eigenvectors are the region’s principal axes, the same computation as principal component analysis of the pixel coordinates, and the angle of the major axis to the $x$ axis is $\theta = \tfrac12 \operatorname{atan2}(2\mu_{11},\, \mu_{20} - \mu_{02})$. A uniformly filled ellipse with semi-axis $a$ has variance $a^2/4$ along it, so the ellipse with the same second moments has full axis lengths $4\sqrt{\lambda_1}$ and $4\sqrt{\lambda_2}$ for the eigenvalues $\lambda_1 \ge \lambda_2$. scikit-image’s regionprops reports exactly these as axis_major_length and axis_minor_length, its eccentricity is $\sqrt{1 - \lambda_2/\lambda_1}$, and its orientation is measured from the row axis rather than the column axis.

Hu’s seven invariants

As Flusser (2000) describes the paper of Hu (1962), Hu introduced moment invariants to the pattern recognition community, used the theory of algebraic invariants to derive seven invariants to rotation, showed how to make them invariant to scaling as well, and demonstrated them on the recognition of printed capital letters. Written with the normalized central moments, as scikit-image and OpenCV compute them,

$$\begin{aligned} \phi_1 &= \eta_{20} + \eta_{02}, \\ \phi_2 &= (\eta_{20} - \eta_{02})^2 + 4\eta_{11}^2, \\ \phi_3 &= (\eta_{30} - 3\eta_{12})^2 + (3\eta_{21} - \eta_{03})^2, \\ \phi_4 &= (\eta_{30} + \eta_{12})^2 + (\eta_{21} + \eta_{03})^2, \\ \phi_5 &= (\eta_{30} - 3\eta_{12})(\eta_{30} + \eta_{12})\big[(\eta_{30} + \eta_{12})^2 - 3(\eta_{21} + \eta_{03})^2\big] \\ &\quad + (3\eta_{21} - \eta_{03})(\eta_{21} + \eta_{03})\big[3(\eta_{30} + \eta_{12})^2 - (\eta_{21} + \eta_{03})^2\big], \\ \phi_6 &= (\eta_{20} - \eta_{02})\big[(\eta_{30} + \eta_{12})^2 - (\eta_{21} + \eta_{03})^2\big] + 4\eta_{11}(\eta_{30} + \eta_{12})(\eta_{21} + \eta_{03}), \\ \phi_7 &= (3\eta_{21} - \eta_{03})(\eta_{30} + \eta_{12})\big[(\eta_{30} + \eta_{12})^2 - 3(\eta_{21} + \eta_{03})^2\big] \\ &\quad - (\eta_{30} - 3\eta_{12})(\eta_{21} + \eta_{03})\big[3(\eta_{30} + \eta_{12})^2 - (\eta_{21} + \eta_{03})^2\big]. \end{aligned}$$

$\phi_1$ and $\phi_2$ use second-order moments only; $\phi_3$ to $\phi_7$ use the third-order ones, $\phi_6$ both. $\phi_1$ is the trace of $\mathbf C$ scaled to unit area, a measure of spread, and $\phi_2$ is the squared difference of its eigenvalues, a measure of elongation. The first six are also unchanged by reflection; $\phi_7$ changes sign, as OpenCV’s documentation of cv::HuMoments states, so it distinguishes a shape from its mirror image.

Complex moments: why the invariants work

The construction is clearest with complex moments, which Flusser (2000) uses to analyse all such invariants:

$$c_{pq} = \iint (x + iy)^p (x - iy)^q f(x, y)\, dx\, dy = \int_0^{2\pi}\!\!\int_0^\infty r^{p+q+1} e^{i(p-q)\theta} f(r, \theta)\, dr\, d\theta .$$

Rotating the image by $\alpha$ multiplies $c_{pq}$ by $e^{-i(p-q)\alpha}$ (Flusser’s Lemma 1, with the rotation written as $f'(r, \theta) = f(r, \theta + \alpha)$): the magnitude stays and only the phase turns. Any product $\prod_j c_{p_j q_j}^{k_j}$ in which the phases cancel, $\sum_j k_j (p_j - q_j) = 0$, is therefore a rotation invariant (his Theorem 1), and with central coordinates and the normalization above it is also invariant to translation and scale. Hu’s invariants are particular products of this kind (Flusser, 2000, eq. 8):

$$\phi_1 = c_{11},\ \phi_2 = c_{20}c_{02},\ \phi_3 = c_{30}c_{03},\ \phi_4 = c_{21}c_{12},\ \phi_5 = \operatorname{Re}(c_{30}c_{12}^3),\ \phi_6 = \operatorname{Re}(c_{20}c_{12}^2),\ \phi_7 = \operatorname{Im}(c_{30}c_{12}^3).$$

Reflection about the $x$ axis replaces every $c_{pq}$ by its complex conjugate, so real parts of such products are unchanged and imaginary parts change sign; Flusser calls the first true invariants and the second skew invariants, the term Hu used (2000, §5). $\phi_7$ is the skew one.

Zernike moments

Zernike polynomials are defined on the unit disk,

$$V_{nl}(\rho, \theta) = R_{nl}(\rho)\, e^{il\theta}, \qquad R_{nl}(\rho) = \sum_{m=0}^{(n-|l|)/2} (-1)^m \frac{(n-m)!}{m!\,\big(\tfrac{n+|l|}{2}-m\big)!\,\big(\tfrac{n-|l|}{2}-m\big)!}\, \rho^{\,n-2m},$$

for $n \ge 0$ and $|l| \le n$ with $n - |l|$ even. They are orthogonal: $\iint_{\rho \le 1} V_{nl} V_{n'l}^{*} = \frac{\pi}{n+1}\delta_{nn'}$. The Zernike moment of order $n$ and repetition $l$ is the projection

$$A_{nl} = \frac{n+1}{\pi} \iint_{\rho \le 1} f(\rho, \theta)\, V_{nl}^{*}(\rho, \theta)\, dx\, dy ,$$

the form that the mahotas library implements with a reference to Teague (1980). Rotation by $\alpha$ shifts only the phase of $A_{nl}$, by $l\alpha$, by the same argument as for complex moments, so the magnitudes $|A_{nl}|$ are rotation invariants. Since $f$ is real, $A_{n,-l} = A_{nl}^{*}$, and orthogonality gives the reconstruction

$$f(x, y) \approx \sum_{n=0}^{N} \sum_{l} A_{nl}\, V_{nl}(x, y).$$

Before the moments are taken, the disk is placed on the shape: its centre at the centroid gives translation invariance, and its radius proportional to a size measure gives scale invariance. The mahotas function leaves the radius to the user and divides by the total mass. The examples below use a radius of $R = 1.2\sqrt{m_{00}}$, our choice, which mirrors the area-based normalization of $\eta_{pq}$ and keeps every pose of the test shape inside the disk.

Flusser (2000) gives the reasons Zernike moments were used: since the polynomials are orthogonal, the moments carry no redundant information and are convenient for reconstructing the image. He reports that Teague (1980) showed the Zernike invariants of second and third order to be equivalent to Hu’s when written in geometric moments, and lists Khotanzad and Hong (1990) among the works that built rotation invariants from Zernike moments.

Properties

Exact in the continuum, approximate on pixels. The invariance proofs treat $f$ as a function of continuous coordinates. A rotated or scaled region on a pixel grid is resampled, its boundary pixels change, and the invariants change with them. OpenCV’s documentation says as much: invariance is proved assuming infinite resolution, and for raster images the computed values of an image and its transformed copy differ a little. The first example below measures how much.

Higher orders are more fragile. A moment of order $p+q$ weights each pixel by distance to the power $p+q$, so the boundary pixels far from the centroid, which are the ones sampling and noise change, dominate higher moments (our reasoning). Flusser (2000, §4.1) advises choosing the fixed moment $c_{p_0 q_0}$ of an invariant basis (below) with the smallest possible indices, because low-order moments are less sensitive to noise than higher-order ones.

Hu’s set is neither independent nor complete. Flusser (2000) showed that a basis of rotation invariants, a smallest set from which all others can be computed, can be built from the products $c_{pq} c_{q_0 p_0}^{p-q}$, $p \ge q$, with one fixed non-zero moment $c_{p_0 q_0}$ with $p_0 - q_0 = 1$, preferably of low order (Theorem 2). Measured against that basis, Hu’s set is dependent: assuming $c_{21} \neq 0$, $\phi_3 = (\phi_5^2 + \phi_7^2)/\phi_4^3$. It is also incomplete: Flusser names two third-order invariants, $c_{20} c_{12}^2$ and $c_{30}^2 c_{02}^3$, that are independent of $\phi_1, \dots, \phi_7$, and his basis of second and third order has six real elements, among them $\operatorname{Im}(c_{20} c_{12}^2)$, which is not one of Hu’s. The fourth-order basis adds five more, eleven in all. In a later note, Flusser (2002) proved that all the moments the basis is built from can be recovered from it (with the rotation fixed by a normalization), so any other rotation invariant of those orders, Hu’s or Zernike’s included, can be computed from it.

Symmetric shapes zero out invariants. If a shape looks the same after a rotation by $2\pi/k$, the rotation property forces $c_{pq} = e^{-i(p-q)2\pi/k}c_{pq}$, so every $c_{pq}$ with $p - q$ not a multiple of $k$ is zero (our derivation from Flusser’s Lemma 1). For a regular five-pointed star all second- and third-order complex moments except $c_{11}$ vanish, and with them $\phi_2$ to $\phi_7$. A shape with a mirror axis equals its own reflection after a rotation, so its skew invariant $\phi_7$ must equal its own negative and is zero. Invariants that are exactly zero carry no information about the shape, and on pixels they come out as small numbers of random size and sign, which the example below shows.

Moments describe the whole region. Every pixel contributes to every moment, so, in our reasoning, a moment changes when any part of the region changes: occlusion, a touching neighbour or a segmentation error far from the centroid shifts all of them. Moment invariants therefore need a clean, isolated region, usually from thresholding and connected components, and they are global shape descriptors rather than local feature descriptors.

Contrast and blur. Moments of a grey-level image depend on its brightness. The normalization of $\eta_{pq}$ is for the size of the region, not its intensity: multiplying the image by a gain $g$ multiplies $\mu_{pq}$ by $g$ and $\mu_{00}^{1+(p+q)/2}$ by $g^{1+(p+q)/2}$, so $\eta_{pq}$ changes by $g^{-(p+q)/2}$, and an offset adds the moments of the whole window (our derivation). Binarizing first removes the dependence on intensity. Flusser (2000) cites moment invariants made invariant also to contrast changes, by Maitra and by Hupkens and de Clippeleir, and moment invariants to linear filtering, by Flusser and colleagues.

Examples

How far the invariants drift on pixels

The script below renders three polygons whose exact moments are known: an L shape with arms of unequal length and no symmetry, an arrow with one mirror axis, and a regular five-pointed star. Each is placed in a $128 \times 128$ image at three sizes (the square root of its area is 10, 20 or 40 pixels), 24 rotations in steps of $15°$ and a random sub-pixel position. Pixel coordinates are $x$ = column and $y$ = row, so $y$ points down; in this left-handed frame $\phi_7$ has the opposite sign from a frame with $y$ up, which affects neither the exact values nor the errors, since both are computed in the same frame.

The exact moments of each posed polygon are computed by Gauss–Legendre quadrature over a fan of triangles, which is exact for the polynomials involved, and give the exact Hu invariants and Zernike moments. The raster ones come from three renderings:

• coverage: each pixel holds the fraction of $8 \times 8$ sub-samples inside the polygon, an anti-aliased grey image;
• binary: a pixel is 1 if its centre is inside, the result of an ideal segmentation;
• noisy (L shape only): the coverage image plus Gaussian noise of standard deviation 0.1, thresholded at 0.5, a segmentation whose boundary pixels flip at random.

For the L shape the script reports the relative error of each Hu invariant and the relative error of the Zernike magnitudes $|A_{nl}|$ of each order $n$ from 2 to 10, taken together over $l$ as $\lVert\, |A_{n\cdot}|_{\text{raster}} - |A_{n\cdot}|_{\text{exact}} \rVert / \lVert\, |A_{n\cdot}|_{\text{exact}} \rVert$. For the arrow and the star it reports the invariants that are exactly zero. Checks come first: the quadrature against a rectangle’s closed-form moments, the exact invariants against rotation and scale, the hand-written Hu invariants against scikit-image’s moments_hu, a finely sampled raster against the exact values, and the properties of the data (which invariants are zero, whether the noisy renders stay one region, whether the shape stays inside the Zernike disk).

The plan written before the code ran stated these expectations: (H1) every median error falls as the shape grows; (H2) the third-order invariants $\phi_3$–$\phi_7$ err more than $\phi_1$ and $\phi_2$; (H3) coverage rendering errs less than binary, and noise more; (H4) the Zernike error rises with the order $n$; (H5) the invariants that are exactly zero come out non-zero, spread over at least a factor of ten across rotations, and the arrow’s $\phi_7$ takes both signs.

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 math import factorial
from scipy import ndimage
from skimage.measure import moments_central, moments_normalized, moments_hu

N_IMG, SIZES, ANGLES, C_DISK = 128, [10, 20, 40], np.arange(24) * 15.0, 1.2
STAR = [((1.0 if k % 2 == 0 else 0.4) * np.cos(np.pi / 2 + k * np.pi / 5),
         (1.0 if k % 2 == 0 else 0.4) * np.sin(np.pi / 2 + k * np.pi / 5)) for k in range(10)]
SHAPES = {"L": np.array([(0, 0), (4, 0), (4, 1), (1, 1), (1, 3), (0, 3)], float),
          "arrow": np.array([(0, -.5), (3, -.5), (3, -1.2), (5, 0), (3, 1.2), (3, .5), (0, .5)], float),
          "star": np.array(STAR, float)}

# --- exact integrals over a polygon: fan of signed triangles, collapsed Gauss-Legendre (exact to degree 78)
gu, gw = np.polynomial.legendre.leggauss(40)
U, W = np.meshgrid((gu + 1) / 2, (gu + 1) / 2, indexing="ij"); WW = np.outer(gw, gw) / 4

def poly_integral(P, g):
    x, y = P[:, 0], P[:, 1]
    sign = np.sign(np.sum(x * np.roll(y, -1) - np.roll(x, -1) * y))
    tot = 0.0
    for b, c in zip(P, np.roll(P, -1, axis=0)):
        s, t = U * (1 - W), U * W                         # triangle (0, b, c)
        tot = tot + (b[0] * c[1] - b[1] * c[0]) * np.sum(WW * U * g(s * b[0] + t * c[0], s * b[1] + t * c[1]))
    return sign * tot

def exact_moments(P):
    m00 = poly_integral(P, lambda x, y: np.ones_like(x))
    xc, yc = poly_integral(P, lambda x, y: x) / m00, poly_integral(P, lambda x, y: y) / m00
    mu = {(p, q): poly_integral(P, lambda x, y, p=p, q=q: (x - xc) ** p * (y - yc) ** q)
          for p in range(4) for q in range(4 - p)}
    return m00, xc, yc, mu

def raster_moments(f):
    yy, xx = np.mgrid[:f.shape[0], :f.shape[1]].astype(float)   # x = column, y = row
    m00 = f.sum(); xc, yc = (f * xx).sum() / m00, (f * yy).sum() / m00
    mu = {(p, q): (f * (xx - xc) ** p * (yy - yc) ** q).sum() for p in range(4) for q in range(4 - p)}
    return m00, xc, yc, mu

def hu(mu):
    e = {k: v / mu[0, 0] ** (1 + sum(k) / 2) for k, v in mu.items()}   # normalized central moments eta_pq
    a, b = e[3, 0] + e[1, 2], e[2, 1] + e[0, 3]
    c, d = e[3, 0] - 3 * e[1, 2], 3 * e[2, 1] - e[0, 3]
    return np.array([e[2, 0] + e[0, 2],
                     (e[2, 0] - e[0, 2]) ** 2 + 4 * e[1, 1] ** 2,
                     c ** 2 + d ** 2,
                     a ** 2 + b ** 2,
                     c * a * (a ** 2 - 3 * b ** 2) + d * b * (3 * a ** 2 - b ** 2),
                     (e[2, 0] - e[0, 2]) * (a ** 2 - b ** 2) + 4 * e[1, 1] * a * b,
                     d * a * (a ** 2 - 3 * b ** 2) - c * b * (3 * a ** 2 - b ** 2)])

# --- Zernike moments (Teague): A_nl = (n+1)/pi * integral over the unit disk of f * conj(V_nl)
def V(n, l, x, y):
    r2, z = x * x + y * y, (x + 1j * y) ** l
    return z * sum((-1) ** m * factorial(n - m) / (factorial(m) * factorial((n + l) // 2 - m)
                   * factorial((n - l) // 2 - m)) * r2 ** ((n - 2 * m - l) // 2) for m in range((n - l) // 2 + 1))

NL = [(n, l) for n in range(2, 11) for l in range(n % 2, n + 1, 2)]

def exact_zernike(P):
    m00, xc, yc, _ = exact_moments(P)
    Q = (P - [xc, yc]) / (C_DISK * np.sqrt(m00))      # unit-disk coordinates
    return {(n, l): (n + 1) / np.pi * poly_integral(Q, lambda x, y: np.conj(V(n, l, x, y))) for n, l in NL}

def raster_zernike(f):
    m00, xc, yc, _ = raster_moments(f); R = C_DISK * np.sqrt(m00)
    yy, xx = np.mgrid[:f.shape[0], :f.shape[1]].astype(float)
    x, y = (xx - xc) / R, (yy - yc) / R; k = x * x + y * y <= 1
    return {(n, l): (n + 1) / np.pi * np.sum(f[k] * np.conj(V(n, l, x[k], y[k]))) / R ** 2 for n, l in NL}

def order_error(Ar, Ae, n):
    vr = np.array([abs(Ar[n, l]) for l in range(n % 2, n + 1, 2)])
    ve = np.array([abs(Ae[n, l]) for l in range(n % 2, n + 1, 2)])
    return np.linalg.norm(vr - ve) / np.linalg.norm(ve)

# --- rendering: k x k sub-samples per pixel (k = 1 is the pixel-centre test)
def render(P, k):
    yy, xx = np.mgrid[:N_IMG, :N_IMG].astype(float); acc = np.zeros((N_IMG, N_IMG))
    for dy in (np.arange(k) + 0.5) / k - 0.5:
        for dx in (np.arange(k) + 0.5) / k - 0.5:
            px, py, inside = xx + dx, yy + dy, np.zeros((N_IMG, N_IMG), bool)
            for (x0, y0), (x1, y1) in zip(P, np.roll(P, -1, axis=0)):   # even-odd crossing test
                cross = (y0 > py) != (y1 > py)
                inside ^= cross & (px < x0 + (py - y0) * (x1 - x0) / np.where(cross, y1 - y0, 1))
            acc += inside
    return acc / k ** 2

def pose(name, size, angle, off):
    P = SHAPES[name]; m00, xc, yc, _ = exact_moments(P)
    a = np.deg2rad(angle); Rm = np.array([[np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]])
    return size / np.sqrt(m00) * (P - [xc, yc]) @ Rm.T + N_IMG / 2 + off

OFFS = np.random.default_rng(0).uniform(-0.5, 0.5, size=(len(SIZES), len(ANGLES), 2))

# ---------------- checks ----------------
m00, xc, yc, mu = exact_moments(np.array([(1, 2), (4, 2), (4, 4), (1, 4)], float))
b1 = max(abs(m00 - 6), abs(xc - 2.5), abs(yc - 3), abs(mu[2, 0] - 4.5), abs(mu[0, 2] - 2), abs(mu[1, 1]), abs(mu[3, 0]))
print(f"B1 rectangle: quadrature vs closed form, max error {b1:.1e} (tol 1e-12)")
EXACT = {}
for si, name in enumerate(SHAPES):
    v = np.array([hu(exact_moments(pose(name, s, a, OFFS[i, j]))[3])
                  for i, s in enumerate(SIZES) for j, a in enumerate(ANGLES)])
    EXACT[name] = v[0]
    if name == "L":
        print(f"B2 L: exact Hu over 72 poses, max relative spread {np.max(np.abs(v - v[0]) / np.abs(v[0])):.1e} (tol 1e-8)")
    else:
        print(f"B2 {name}: exact Hu over 72 poses, max spread / max|phi| {np.max(np.abs(v - v[0])) / np.abs(v[0]).max():.1e} (tol 1e-10)")
fb = render(pose("L", 20, 30.0, OFFS[1, 2]), 1)
sk = moments_hu(moments_normalized(moments_central(fb.T, order=3), order=3))   # transposed: axis 0 is x
print(f"B3 hand Hu vs skimage moments_hu, max relative difference {np.max(np.abs(hu(raster_moments(fb)[3]) - sk) / np.abs(sk)):.1e} (tol 1e-10)")
A0 = None; b5 = 0.0
for i, s in enumerate(SIZES):
    for j in range(0, 24, 6):
        A = exact_zernike(pose("L", s, ANGLES[j], OFFS[i, j])); mag = np.array([abs(A[k]) for k in NL])
        A0 = mag if A0 is None else A0; b5 = max(b5, np.max(np.abs(mag - A0)) / A0.max())
print(f"B5 exact |A_nl| (n<=10) over 12 poses, max spread / max|A| {b5:.1e} (tol 1e-10)")
P6 = pose("L", 40, 30.0, OFFS[2, 2]); f6 = render(P6, 32)
r6, e6 = raster_moments(f6), exact_moments(P6)
print(f"B6 fine raster (32x32 sub-samples, size 40) vs exact: m00 {abs(r6[0] - e6[0]) / e6[0]:.1e} (tol 1e-5), "
      f"Hu max {np.max(np.abs(hu(r6[3]) - hu(e6[3])) / np.abs(hu(e6[3]))):.1e} (tol 1e-3), Zernike orders 2/6/10 "
      + "/".join(f"{order_error(raster_zernike(f6), exact_zernike(P6), n):.1e}" for n in (2, 6, 10)) + " (tol 1e-2)")
print(f"D1 L: min|phi|/max|phi| = {np.abs(EXACT['L']).min() / np.abs(EXACT['L']).max():.1e} (>= 1e-8: no exact zero)")
print(f"D2 arrow: |phi7|/max|phi| = {abs(EXACT['arrow'][6]) / np.abs(EXACT['arrow']).max():.1e} (<= 1e-12)")
print(f"D3 star: max|phi2..7|/phi1 = {np.abs(EXACT['star'][1:]).max() / EXACT['star'][0]:.1e} (<= 1e-12)")

# ---------------- example: raster invariants against the exact ones ----------------
CONDS = ["coverage", "binary", "noisy"]
herr = {(c, s): [] for c in CONDS for s in SIZES}; zerr = {(c, s): [] for c in CONDS for s in SIZES}
zero_vals = {(nm, s): [] for nm in ("arrow", "star") for s in SIZES}
multi = outside = 0; dmax = 0.0
for i, s in enumerate(SIZES):
    for j, a in enumerate(ANGLES):
        P = pose("L", s, a, OFFS[i, j]); m00, xc, yc, mu = exact_moments(P)
        dmax = max(dmax, np.max(np.hypot(*(P - [xc, yc]).T)) / (C_DISK * np.sqrt(m00)))
        h_ex, z_ex = hu(mu), exact_zernike(P)
        fc = render(P, 8)
        fn = (fc + np.random.default_rng(100 * i + j).normal(0, 0.1, fc.shape) > 0.5).astype(float)
        multi += ndimage.label(fn)[1] != 1
        outside += np.sum((fn > 0) & ~ndimage.binary_dilation(fc > 0))
        for c, f in zip(CONDS, (fc, render(P, 1), fn)):
            herr[c, s].append(np.abs(hu(raster_moments(f)[3]) - h_ex) / np.abs(h_ex))
            zr = raster_zernike(f); zerr[c, s].append([order_error(zr, z_ex, n) for n in range(2, 11)])
        for nm in ("arrow", "star"):
            zero_vals[nm, s].append(hu(raster_moments(render(pose(nm, s, a, OFFS[i, j]), 1))[3]))
print(f"D4 L noisy renders: {multi} with != 1 component, {outside} foreground pixels off the true support (both 0)")
print(f"D5 L: max vertex distance / Zernike disk radius = {dmax:.3f} (< 1)")

print("\nExact Hu invariants (phi1..phi7):")
for nm, v in EXACT.items():
    print(f"  {nm:5s} " + " ".join(f"{x: .2e}" for x in v))
hdr = "".join(f"{c[:3]}{s:>3d}" .rjust(8) for c in CONDS for s in SIZES)
for stat, fn_ in (("median", np.median), ("max", np.max)):
    print(f"\nL shape, Hu relative error (%), {stat} over 24 rotations; columns: condition, sqrt(area) px")
    print("      " + hdr)
    for k in range(7):
        print(f"  phi{k + 1}" + "".join(f"{100 * fn_(np.array(herr[c, s])[:, k]):8.3f}" for c in CONDS for s in SIZES))
print("\nL shape, Zernike order-n magnitude error (%), median over 24 rotations")
print("      " + hdr)
for n in range(2, 11):
    print(f"  n={n:<3d}" + "".join(f"{100 * np.median(np.array(zerr[c, s])[:, n - 2]):8.3f}" for c in CONDS for s in SIZES))

print("\nInvariants that are exactly zero, binary renders: log10|phi| min..max over 24 rotations")
print("  (L's exact log10|phi| for scale: " + " ".join(f"{np.log10(abs(x)):.1f}" for x in EXACT["L"]) + ")")
print("              " + "".join(f"size {s:<2d}".rjust(14) for s in SIZES))
rows = [("arrow", 6)] + [("star", k) for k in range(1, 7)]
for nm, k in rows:
    cells = []
    for s in SIZES:
        lv = np.log10(np.abs(np.array(zero_vals[nm, s])[:, k]))
        cells.append(f"{lv.min():6.1f}..{lv.max():5.1f}")
    print(f"  {nm:5s} phi{k + 1}  " + "".join(c.rjust(14) for c in cells))
print("  arrow phi7 signs (+/-) per size: "
      + ", ".join(f"{np.sum(np.array(zero_vals['arrow', s])[:, 6] > 0)}/{np.sum(np.array(zero_vals['arrow', s])[:, 6] < 0)}" for s in SIZES))

B1 rectangle: quadrature vs closed form, max error 8.0e-15 (tol 1e-12)
B2 L: exact Hu over 72 poses, max relative spread 1.5e-10 (tol 1e-8)
B2 arrow: exact Hu over 72 poses, max spread / max|phi| 3.9e-13 (tol 1e-10)
B2 star: exact Hu over 72 poses, max spread / max|phi| 1.0e-12 (tol 1e-10)
B3 hand Hu vs skimage moments_hu, max relative difference 8.9e-13 (tol 1e-10)
B5 exact |A_nl| (n<=10) over 12 poses, max spread / max|A| 1.7e-13 (tol 1e-10)
B6 fine raster (32x32 sub-samples, size 40) vs exact: m00 6.1e-07 (tol 1e-5), Hu max 3.0e-04 (tol 1e-3), Zernike orders 2/6/10 2.7e-04/8.9e-04/3.3e-03 (tol 1e-2)
D1 L: min|phi|/max|phi| = 1.6e-05 (>= 1e-8: no exact zero)
D2 arrow: |phi7|/max|phi| = 6.6e-18 (<= 1e-12)
D3 star: max|phi2..7|/phi1 = 2.2e-25 (<= 1e-12)
D4 L noisy renders: 0 with != 1 component, 0 foreground pixels off the true support (both 0)
D5 L: max vertex distance / Zernike disk radius = 0.916 (< 1)

Exact Hu invariants (phi1..phi7):
  L      3.47e-01  4.34e-02  2.31e-02  2.57e-03 -5.56e-06 -1.50e-04 -1.91e-05
  arrow  3.38e-01  7.94e-02  3.18e-03  1.50e-03  3.28e-06  4.23e-04 -2.22e-18
  star   2.10e-01  9.83e-28  4.67e-26  1.41e-26  2.85e-52 -4.09e-40  2.19e-52

L shape, Hu relative error (%), median over 24 rotations; columns: condition, sqrt(area) px
        cov 10  cov 20  cov 40  bin 10  bin 20  bin 40  noi 10  noi 20  noi 40
  phi1   0.477   0.121   0.031   0.948   0.280   0.077   1.873   0.645   0.172
  phi2   0.071   0.020   0.008   3.900   1.072   0.297   4.358   2.360   0.639
  phi3   0.088   0.029   0.005   3.153   0.961   0.284   8.798   2.566   0.790
  phi4   0.219   0.034   0.009   4.540   1.331   0.417   8.299   2.602   1.241
  phi5   1.104   0.173   0.042  29.343   8.742   1.914  28.741  21.020   5.003
  phi6   0.960   0.150   0.041  24.382   6.666   1.415  30.512  19.101   3.405
  phi7   0.312   0.054   0.017   8.292   2.311   0.848  17.954   6.365   2.546

L shape, Hu relative error (%), max over 24 rotations; columns: condition, sqrt(area) px
        cov 10  cov 20  cov 40  bin 10  bin 20  bin 40  noi 10  noi 20  noi 40
  phi1   1.059   0.482   0.434   5.001   4.209   1.965   5.365   3.963   2.463
  phi2   2.627   1.193   0.985  22.084  12.953   3.914  11.757  11.381   4.775
  phi3   3.893   2.315   1.756  29.450  17.860   9.444  23.883  15.776  10.725
  phi4   5.014   3.579   1.810  64.330  27.766  15.536  35.354  23.448  13.731
  phi5  16.248   8.839   5.564 256.865  96.547  63.319 131.627  71.105  48.629
  phi6  11.632   6.870   4.317 119.998  63.995  43.881  91.183  51.147  33.765
  phi7  10.121   7.359   3.458 145.436  54.256  31.476  76.099  48.537  26.163

L shape, Zernike order-n magnitude error (%), median over 24 rotations
        cov 10  cov 20  cov 40  bin 10  bin 20  bin 40  noi 10  noi 20  noi 40
  n=2     0.431   0.108   0.028   1.051   0.286   0.082   1.775   0.632   0.179
  n=3     0.085   0.016   0.004   2.415   0.554   0.181   4.418   1.352   0.538
  n=4     0.949   0.241   0.060   3.285   1.431   0.248   4.619   2.296   0.696
  n=5     1.458   0.363   0.093   2.082   0.640   0.195   2.910   1.159   0.358
  n=6     1.446   0.374   0.091   7.399   1.952   0.419   9.487   3.237   1.042
  n=7     2.310   0.579   0.145   3.432   1.022   0.241   4.832   1.856   0.536
  n=8     2.445   0.605   0.149   7.079   2.898   0.769  11.050   4.222   1.444
  n=9     3.638   0.936   0.235   7.609   1.853   0.413   8.122   2.788   0.898
  n=10    5.029   1.286   0.328  15.304   3.161   0.874  15.012   5.652   1.711

Invariants that are exactly zero, binary renders: log10|phi| min..max over 24 rotations
  (L's exact log10|phi| for scale: -0.5 -1.4 -1.6 -2.6 -5.3 -3.8 -4.7)
                     size 10       size 20       size 40
  arrow phi7    -23.7.. -6.4   -9.1.. -6.8   -inf.. -7.3
  star  phi2     -5.9.. -3.7   -6.9.. -4.4   -7.6.. -5.2
  star  phi3     -5.7.. -3.8   -8.2.. -4.4   -7.9.. -5.3
  star  phi4     -6.6.. -4.4   -7.7.. -5.6   -7.9.. -6.2
  star  phi5    -13.4.. -9.1  -14.5..-11.0  -16.8..-12.3
  star  phi6     -9.6.. -6.7  -11.3.. -8.0  -12.4.. -9.2
  star  phi7    -13.2.. -9.3  -15.0..-10.9  -15.5..-12.7
  arrow phi7 signs (+/-) per size: 14/10, 11/13, 11/12

All checks pass, so the differences in the tables are due to the raster alone: the exact invariants agree across the 72 poses to rounding error, the hand-written invariants match scikit-image, and a render with $32 \times 32$ sub-samples per pixel matches the exact values to within 0.03% for the Hu invariants and 0.33% for the Zernike magnitudes of orders 2, 6 and 10.

Size (H1 holds). Every median error falls as the shape grows, under all three conditions, in both tables. From one size to the next it falls by a factor of between about 2.3 and 6.4 under coverage and binary rendering (our arithmetic on the printed medians), so a shape twice as wide is measured several times more accurately. At size 10, the binary $\phi_5$ and $\phi_6$ are off by a median of 24–29%.

Order (H2 holds). In every column the largest error among $\phi_3$–$\phi_7$ exceeds that of $\phi_1$ and $\phi_2$. The worst are $\phi_5$ and $\phi_6$, with medians up to 29% and maxima up to 257% for binary renders at size 10. These two are also among the smallest invariants of the L shape ($10^{-5.3}$ and $10^{-3.8}$, against $10^{-0.5}$ for $\phi_1$; $\phi_7$ is $10^{-4.7}$). In our reading, each is a difference of products of third-order moments that nearly cancel for this shape, so a small absolute error in the moments becomes a large relative error; an invariant’s reliability depends on the shape as well as its order.

Rendering and noise (H3: the coverage part holds, the noise part fails). Coverage rendering errs less than binary for every invariant and size, in both tables, by a factor of about two to more than fifty: the grey boundary pixels record how much of each pixel the shape covers, while the binary image moves the boundary by up to half a pixel. Noise increases the error over binary in every entry but one: at size 10 the median $\phi_5$ error is 28.7% with noise against 29.3% without, and in the Zernike table the order-10 error at size 10 is 15.0% with noise against 15.3% without. So the noise part of H3 fails. In both exceptions the errors are large, and the other 20 Hu medians and 26 Zernike medians rise with noise. The maxima over rotations are many times the medians: a few poses are much worse than the others, so a shape classifier tuned on a few rotations can be surprised by a new one.

Zernike order (H4 fails). Zernike errors grow with order overall, from 0.43% at $n = 2$ to 5.0% at $n = 10$ for coverage at size 10, but not at every step as planned. Under binary and noisy rendering the odd orders 5, 7 and, in five of the six columns, 9 come out below the even order before them (at size 10, binary: 3.3% at $n = 4$, 2.1% at $n = 5$, 7.4% at $n = 6$); under coverage the error drops at $n = 3$ and, at sizes 10 and 40, at $n = 6$. The run does not show why. For the binary and noisy columns, one untested possibility, our conjecture: only even orders contain the rotationally symmetric term $l = 0$, which weighs the whole region’s radial mass profile and may be the most affected by the radius set from the raster area. At the second order the Zernike error is close to that of $\phi_1$, consistent, in our reading, with the equivalence of the second- and third-order invariants of the two families that Flusser (2000) reports from Teague (1980).

Symmetric shapes (H5 holds except at one pose). Invariants that are exactly zero come out non-zero on every raster but one, at magnitudes that range over about 1.7 powers of ten or more across the 24 rotations. The star’s $\phi_2$ reaches $10^{-3.7}$ at size 10, small compared with the L shape’s $10^{-1.4}$ but not zero, and its logarithm, the form in which OpenCV’s cv::matchShapes compares Hu invariants, varies between $-5.9$ and $-3.7$ for one shape at one size. The arrow’s $\phi_7$ takes both signs at every size (14 positive and 10 negative at size 10). At size 40 one rotation gives exactly zero, printed as $-\infty$, which is why the signs there add up to 23, and H5’s expectation of a non-zero value fails for that pose; a log transform of its $\phi_7$ is undefined. For symmetric shapes, in our reading, the vanishing invariants are better left out or replaced by invariants of orders that the symmetry does not cancel.

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 3 seconds on one CPU, and two runs printed identical output.

Reconstruction from Zernike moments

Orthogonality means that the Zernike moments of orders up to $N$ give the best approximation of the image by Zernike polynomials of those orders, in the least-squares sense over the disk, so they can be checked by reconstructing the image. The second script renders the L shape with coverage at a size that gives a disk radius of 56 pixels, computes $A_{nl}$ up to order 32 in two ways, by pixel sums and by the exact polygon integrals, and reconstructs the image at pixel centres for $N$ from 2 to 32. It reports the relative $L_2$ error over the disk and the number of misclassified pixels after thresholding at 0.5, as a fraction of the shape’s pixels.

It first checks the polynomials. The direct factorial formula for $R_{nl}$ was found during planning to lose accuracy above order 32, which is why $N$ stops there; the script prints the orthogonality error at orders up to 40 and 48 for information. It also builds an image from six known Zernike coefficients, sampled on the same disk, and recovers them by pixel sums.

The plan’s expectations: (H6) with exact coefficients the error decreases steadily as $N$ grows; (H7) with coefficients from pixel sums the error is close to the exact-coefficient error at low $N$ and larger at the highest $N$.

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 math import factorial

N_IMG, C_DISK = 128, 1.2
L_SHAPE = np.array([(0, 0), (4, 0), (4, 1), (1, 1), (1, 3), (0, 3)], float)

gu, gw = np.polynomial.legendre.leggauss(40)
U, W = np.meshgrid((gu + 1) / 2, (gu + 1) / 2, indexing="ij"); WW = np.outer(gw, gw) / 4

def poly_integral(P, g):                      # exact for polynomials of degree <= 78
    x, y = P[:, 0], P[:, 1]
    sign = np.sign(np.sum(x * np.roll(y, -1) - np.roll(x, -1) * y))
    tot = 0.0
    for b, c in zip(P, np.roll(P, -1, axis=0)):
        s, t = U * (1 - W), U * W
        tot = tot + (b[0] * c[1] - b[1] * c[0]) * np.sum(WW * U * g(s * b[0] + t * c[0], s * b[1] + t * c[1]))
    return sign * tot

def V(n, l, x, y):                            # Zernike polynomial, direct factorial formula
    r2, z = x * x + y * y, (x + 1j * y) ** l
    return z * sum((-1) ** m * factorial(n - m) / (factorial(m) * factorial((n + l) // 2 - m)
                   * factorial((n - l) // 2 - m)) * r2 ** ((n - 2 * m - l) // 2) for m in range((n - l) // 2 + 1))

def pairs(nmax):
    return [(n, l) for n in range(nmax + 1) for l in range(n % 2, n + 1, 2)]

def render(P, k):
    yy, xx = np.mgrid[:N_IMG, :N_IMG].astype(float); acc = np.zeros((N_IMG, N_IMG))
    for dy in (np.arange(k) + 0.5) / k - 0.5:
        for dx in (np.arange(k) + 0.5) / k - 0.5:
            px, py, inside = xx + dx, yy + dy, np.zeros((N_IMG, N_IMG), bool)
            for (x0, y0), (x1, y1) in zip(P, np.roll(P, -1, axis=0)):
                cross = (y0 > py) != (y1 > py)
                inside ^= cross & (px < x0 + (py - y0) * (x1 - x0) / np.where(cross, y1 - y0, 1))
            acc += inside
    return acc / k ** 2

# ---------------- checks ----------------
gr, gwr = np.polynomial.legendre.leggauss(80)
rq, wr = (gr + 1) / 2, gwr / 2; tq = 2 * np.pi * np.arange(160) / 160
RQ, TQ = np.meshgrid(rq, tq, indexing="ij"); WQ = np.outer(wr, np.full(160, 2 * np.pi / 160)) * RQ
XQ, YQ = RQ * np.cos(TQ), RQ * np.sin(TQ)
for nmax in (16, 32, 40, 48):
    Vq = {p: V(*p, XQ, YQ) for p in pairs(nmax)}; worst = 0.0
    for (n1, l1) in Vq:
        for (n2, l2) in Vq:
            if l1 == l2 and n2 >= n1:
                ip = np.sum(WQ * Vq[n1, l1] * np.conj(Vq[n2, l2]))
                worst = max(worst, abs(ip - (np.pi / (n1 + 1) if n1 == n2 else 0)) * (n1 + 1) / np.pi)
    tol = {16: " (tol 1e-10)", 32: " (tol 1e-5)"}.get(nmax, " (information only)")
    print(f"B4 orthogonality of V_nl by quadrature, n <= {nmax}: max error {worst:.1e}{tol}")

# the shape: L, sqrt(area) = 56/1.2 px, rotated 30 degrees, centroid at pixel (64, 64)
a = np.deg2rad(30.0); Rm = np.array([[np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]])
m00 = poly_integral(L_SHAPE, lambda x, y: np.ones_like(x))
cen = [poly_integral(L_SHAPE, lambda x, y: x) / m00, poly_integral(L_SHAPE, lambda x, y: y) / m00]
P = (56 / C_DISK) / np.sqrt(m00) * (L_SHAPE - cen) @ Rm.T + N_IMG / 2
f = render(P, 8)
yy, xx = np.mgrid[:N_IMG, :N_IMG].astype(float)
xc, yc = (f * xx).sum() / f.sum(), (f * yy).sum() / f.sum(); R = C_DISK * np.sqrt(f.sum())
x, y = (xx - xc) / R, (yy - yc) / R; k = x * x + y * y <= 1; x, y, fk = x[k], y[k], f[k]
print(f"disk: centre ({xc:.3f}, {yc:.3f}), R = {R:.3f} px, {k.sum()} pixels; "
      f"max vertex distance / R = {np.max(np.hypot(P[:, 0] - xc, P[:, 1] - yc)) / R:.3f}")

known = {(0, 0): .5, (2, 0): .2, (2, 2): .1 + .05j, (3, 1): -.08j, (4, 4): .06, (6, 2): .04 - .03j}
g = sum((c * V(n, l, x, y) * (1 if l == 0 else 2)).real for (n, l), c in known.items())
est = {p: (p[0] + 1) / np.pi * np.sum(g * np.conj(V(*p, x, y))) / R ** 2 for p in pairs(10)}
print(f"B7 known coefficients recovered by pixel sums: max error {max(abs(est[p] - c) for p, c in known.items()):.1e} "
      f"(tol 1.5e-2), max leak into absent terms {max(abs(v) for p, v in est.items() if p not in known):.1e} (tol 3e-2)")

# ---------------- reconstruction from raster and from exact coefficients ----------------
Q = (P - [xc, yc]) / R
Vk = {p: V(*p, x, y) for p in pairs(32)}
A_r = {p: (p[0] + 1) / np.pi * np.sum(fk * np.conj(v)) / R ** 2 for p, v in Vk.items()}
A_e = {p: (p[0] + 1) / np.pi * poly_integral(Q, lambda u, w, p=p: np.conj(V(*p, u, w))) for p in Vk}
inside = fk >= 0.5
print("\n  N  terms  L2 error (raster A)  L2 error (exact A)  misclassified/shape px (raster A, exact A)")
for N in (2, 4, 8, 12, 16, 20, 24, 28, 32):
    out = []
    for A in (A_r, A_e):
        fh = sum((A[n, l] * Vk[n, l] * (1 if l == 0 else 2)).real for n, l in pairs(N))
        out += [np.linalg.norm(fh - fk) / np.linalg.norm(fk), np.sum((fh > 0.5) != inside) / inside.sum()]
    print(f"{N:3d} {len(pairs(N)):6d} {out[0]:16.4f} {out[2]:19.4f} {out[1]:16.4f} {out[3]:9.4f}")

B4 orthogonality of V_nl by quadrature, n <= 16: max error 8.9e-13 (tol 1e-10)
B4 orthogonality of V_nl by quadrature, n <= 32: max error 7.5e-07 (tol 1e-5)
B4 orthogonality of V_nl by quadrature, n <= 40: max error 7.3e-04 (information only)
B4 orthogonality of V_nl by quadrature, n <= 48: max error 7.4e+01 (information only)
disk: centre (64.000, 64.000), R = 56.000 px, 9844 pixels; max vertex distance / R = 0.916
B7 known coefficients recovered by pixel sums: max error 7.3e-03 (tol 1.5e-2), max leak into absent terms 1.5e-02 (tol 3e-2)

  N  terms  L2 error (raster A)  L2 error (exact A)  misclassified/shape px (raster A, exact A)
  2      4           0.7513              0.7513           0.8867    0.8867
  4      9           0.6691              0.6691           0.5876    0.5876
  8     25           0.4588              0.4588           0.1784    0.1784
 12     49           0.3308              0.3308           0.0674    0.0674
 16     81           0.3010              0.3010           0.0610    0.0610
 20    121           0.2848              0.2848           0.0427    0.0417
 24    169           0.2528              0.2528           0.0317    0.0321
 28    225           0.2169              0.2169           0.0234    0.0239
 32    289           0.1966              0.1967           0.0193    0.0197

The polynomials are orthogonal to $10^{-12}$ up to order 16 and to $10^{-6}$ up to order 32. Above that the direct formula breaks down, with an error of $7 \times 10^{-4}$ up to order 40 and 74 up to order 48: its terms alternate in sign and grow with the factorials, so in our reading floating-point cancellation destroys the result, and higher orders need a numerically stable way to evaluate $R_{nl}$, which this code does not implement. On a disk of 56 pixels, pixel sums recover known coefficients to within 0.007, with leakage of up to 0.015 into absent terms: the sampled polynomials are close to orthogonal, not exactly.

H6 holds. With exact coefficients the $L_2$ error falls at every step, from 0.75 at $N = 2$ to 0.20 at $N = 32$, and the misclassified fraction falls from 89% to 2%. At low orders the reconstruction is a smooth blob that cannot be thresholded into an L; by order 12 the threshold recovers the shape to within 7% of its pixels. The $L_2$ error falls slowly because the shape’s straight edges and corners are discontinuities, which a finite sum of smooth polynomials approximates with overshoot and ringing near the edges, the same effect as truncating a Fourier series (our reading).

H7 fails. Pixel-sum coefficients give the same $L_2$ errors as the exact coefficients to four decimals up to $N = 28$ and the same misclassified fractions up to $N = 16$; beyond that they differ by at most 0.001: at $N = 32$ the pixel-sum reconstruction is in fact very slightly better ($L_2$ error 0.1966 against 0.1967, 1.93% against 1.97% misclassified). On a disk of this size, up to order 32, the limit on reconstruction is the truncation of the series, not the sampling of the moments. In our reading, the pixel sums may do slightly better because they are computed from the same pixel grid on which the error is measured, while the exact coefficients describe the continuous polygon.

The script ran in about 1 second. Both examples ran as piloted, with no change after the first run. Two planning decisions made before the pilot shaped them: the cap at order 32, set from the orthogonality check shown above, and applying the noisy condition to the L shape only, for which check D4 confirms that every noisy render stays one region.

Common Misconceptions

• “Hu moments are invariant, so the same shape always gives the same values.” They are invariant for continuous shapes. On pixels the third-order invariants of a binary shape 10 pixels across erred by a median of up to 29% over rotations in the example, and by more than 100% for some poses. One remedy, comparing the logarithms of the invariants, as OpenCV’s cv::matchShapes does with $\operatorname{sign}(h_i)\log_{10}|h_i|$, compresses the range but does not remove the error.
• “Seven Hu invariants describe a shape completely.” They are dependent ($\phi_3$ follows from $\phi_4$, $\phi_5$ and $\phi_7$) and incomplete (Flusser, 2000), and they use only moments up to third order, which cannot tell apart shapes that differ in finer details. Different shapes can have the same seven values.
• “A small invariant is a reliable feature.” For symmetric shapes, invariants that are zero in theory are pure sampling noise on pixels, and their logarithms swing by orders of magnitude.
• “Moments are a contour feature.” Moments of a filled region and of its outline are different quantities. OpenCV’s cv::moments computes the moments of the polygon when given contour points, by a sum over its edges, and those of the pixel values when given an image; scikit-image computes moments of an image (moments) or of a set of points (moments_coords).
• “Normalized moments make features independent of brightness.” $\eta_{pq}$ is normalized for the region’s size: on a grey-level image a gain $g$ scales it by $g^{-(p+q)/2}$, and an offset or a non-uniform background adds moments of its own (Properties). Binarizing first removes the dependence on intensity but adds the thresholding error shown in the example.
• “Higher-order Zernike moments are always better for reconstruction.” In the example, the error fell steadily up to order 32, but the direct formula for the polynomials breaks down above that, and the edge ringing persists at every finite order.

Where It Is Used

• Region properties. Area, centroid, orientation, axis lengths and eccentricity of segmented regions, as in scikit-image’s regionprops, are moments of order up to two; they are measured after image segmentation in microscopy and industrial inspection.
• Tracking. CAMShift sets its window’s size and orientation from the moments of a back-projected probability image; Mean-Shift Tracking covers it.
• Keypoint orientation. ORB orients each corner by the intensity centroid of its patch, the direction from the keypoint to $(m_{10}/m_{00}, m_{01}/m_{00})$.
• Oriented boxes and pose. Second moments fit an oriented box to a segmented object in aerial object detection and give a rough position, scale and angle to start a fit, as in one of the facial landmark detection baselines.
• Visual servoing. Image moments of planar objects of any shape serve as visual features with known interaction matrices in visual servoing.
• Shape classification. Flusser (2000) calls moment invariants a classical tool for object recognition; in the semi-supervised learning example, the logarithms of the first three Hu invariants describe shapes for a linear classifier. Flusser (2000) lists early applications to aircraft silhouettes, character recognition, registration of satellite images and industrial inspection.
• Normalization. Moments define a canonical position, scale and orientation; Flusser (2000) cites their use for the geometric normalization of images. The deslanted variant of MNIST shears each digit so that the principal axis of its ink, computed from its second moments, is vertical.