An Edge-Based Formulation for the Exact Computation of High-Order Zernike Moments of 2D Shapes and Images

Abstract: Zernike moments are widely used rotation-invariant descriptors for shape and image analysis, but their standard computation relies on a pixel-based quadrature that treats each pixel as a point mass located at its center. This approximation introduces spatial aliasing that increases with moment order, degrading image reconstruction and reducing the discriminative power of high-order moments. We present an edge-based formulation that eliminates this source of error by applying Green's theorem to transform the two-dimensional area integral defining a Zernike moment into a sum of one-dimensional integrals along image boundaries. The resulting framework applies equally to polygonal shapes, binary images, grayscale images, and color images. We derive recurrence relations for the required radial primitives and show that the transformed edge integrands are polynomial functions, allowing their exact evaluation using Clenshaw--Curtis quadrature. The proposed method computes Zernike moments from polygonal image representations without the spatial discretization errors inherent to conventional pixel-based approaches and remains computationally practical for high-order moments. Numerical experiments on image reconstruction, shape analysis, and character classification demonstrate that the proposed formulation matches the accuracy of classical methods at low orders while remaining stable at orders for which pixel-based moments suffer from significant aliasing and numerical degradation.

• The paper transforms 2D Zernike moment area integrals into boundary line integrals and proves that edge integrands are polynomials of degree at most n, enabling exact Clenshaw–Curtis evaluation with at least n + 1 nodes.
• The method supports polygonal contours, binary, grayscale, raster, and quaternion color images while reducing active computation to image edges and achieving practical near-linear scaling with the number of moments.
• Experiments show that edge-based moments prevent high-order aliasing, preserve QR-code reconstruction, and maintain 97.3% classification accuracy at order 500 versus 90.1% for pixel-based moments, though computation is 37–74 times slower.

Overview and motivation

Koehl and Tillmann address a long-standing numerical deficiency in the computation of 2D Zernike moments: the conventional pixel-based quadrature, which models each pixel as a Dirac point mass located at its center. This lowest-order quadrature introduces spatial aliasing that grows with moment order, since the radial polynomials $R_n^{|m|}(r)$ oscillate on a scale comparable to a single pixel once the order is sufficiently high. The authors propose an alternative in which the two-dimensional area integral defining a Zernike moment is transformed, via Green's theorem, into a sum of one-dimensional line integrals along image boundaries. The resulting formulation is exact for the piecewise-constant image model, applies uniformly to polygonal contours, binary images, grayscale images, and color images, and eliminates discretization error at its source rather than correcting it after the fact.

Mathematical formulation

The method rests on a decomposition of the domain $D \subseteq \mathbb{D}$ (the unit disk) into signed triangles $\sigma_{AB} = (O, A, B)$ formed by the origin $O$ and each directed edge $e = (A, B)$ of the boundary, with orientation factors determined by the signed area $\tfrac{1}{2}\det(\mathbf{A}, \mathbf{B})$. For each edge, the area integral over the associated triangle is separated in polar coordinates into a radial primitive

$Q_n^m(r) = \int_0^r t\, R_n^m(t)\, \mathrm{d}t$

followed by an angular line integral $$Z_n^m(e) = \epsilon_{AB}\int_{\theta_A}^{\theta_B} Q_n^m(r(\theta))\, e^{-im\theta}\,\mathrm{d}\theta$$. Recurrences for $Q_n^m$ and the auxiliary primitive $S_n^m(r) = \int_0^r R_n^m(t)\,\mathrm{d}t$ are derived from Bessel-function integral representations of the Zernike radial polynomials; the recurrence for $D \subseteq \mathbb{D}$0 is identified as new, while the others restate known results with complete proofs supplied in the appendix.

The central theoretical result is that, after parameterizing the edge linearly by $D \subseteq \mathbb{D}$1, the transformed integrand

$D \subseteq \mathbb{D}$2

is a polynomial in $D \subseteq \mathbb{D}$3 of degree at most $D \subseteq \mathbb{D}$4. The proof exploits the Jacobi-polynomial representation $D \subseteq \mathbb{D}$5, which yields $D \subseteq \mathbb{D}$6, and the identity $D \subseteq \mathbb{D}$7. Because the integrand is a polynomial of known degree, the edge integral can be evaluated exactly using Clenshaw–Curtis quadrature with $D \subseteq \mathbb{D}$8 Chebyshev–Gauss–Lobatto nodes. The choice of Clenshaw–Curtis over Gauss–Legendre is justified by the nested structure of its nodes, which permits adaptive refinement without discarding prior evaluations and allows a single set of node evaluations to be shared across all moments up to the maximum order.

Generalization and algorithm

The framework extends beyond closed contours to arbitrary 2D meshes: for a polygonal tessellation, contributions along internal edges cancel by orientation, so only true external boundaries and hole boundaries contribute. Applied to a raster image treated as a uniform square mesh, each edge's net contribution is weighted by the directional intensity gradient $D \subseteq \mathbb{D}$9 across it; edges between identically valued pixels vanish. This gradient weighting implies that for structured images the active edge count scales with the shape perimeter rather than the image area.

Color images are handled through right-sided Quaternion Zernike Moments, with pixels encoded as pure quaternions $\sigma_{AB} = (O, A, B)$0 and the complex exponential replaced by a quaternion exponential about the space-diagonal axis $\sigma_{AB} = (O, A, B)$1. The authors note a practical asymmetry here: unlike the complex case, where $\sigma_{AB} = (O, A, B)$2 halves the work, quaternion moments admit no such conjugation shortcut because quaternion multiplication is non-commutative.

On complexity, the exact edge-based computation is $\sigma_{AB} = (O, A, B)$3 in the worst case, compared with $\sigma_{AB} = (O, A, B)$4 for an exact area-based pixel quadrature and $\sigma_{AB} = (O, A, B)$5 for the (inexact) point-mass method. In practice, adaptive Clenshaw–Curtis quadrature converged with an average of 17 nodes per edge on $\sigma_{AB} = (O, A, B)$6 images, making observed runtime nearly linear in the number of moments $\sigma_{AB} = (O, A, B)$7. The edge loop is embarrassingly parallel.

Experimental findings

Three classes of experiments compare the edge-based method against the standard pixel-based approach, using C++ implementations (Shape2Zernike and Zernike2Shape) on a 16-core Apple M4 Max.

Reconstruction fidelity. For the $\sigma_{AB} = (O, A, B)$8 Barbara color image and a 64-tooth grayscale gear image, quaternion and complex moments up to order $\sigma_{AB} = (O, A, B)$9 were computed with both methods. Reconstructions at orders $O$0 and $O$1 were qualitatively and quantitatively indistinguishable between methods, with contour-based Hausdorff distances of one pixel at $O$2 for all three computation variants. This establishes that the proposed method preserves the accuracy of the classical approach at moderate orders on high-resolution imagery.

The cost of exactness. The authors report that edge-based computation is 37 times slower than pixel-based for the gear image and 74 times slower for Barbara, attributable to roughly two unique edges per pixel multiplied by ~17 quadrature evaluations each. They state plainly that this one-to-two order-of-magnitude slowdown motivates the question of whether the extra accuracy matters — and answer it with the aliasing experiments below.

Ghost-pixel artifacts. On a $O$3 binary image containing a single white pixel, pixel-based moments produce a spurious "ghost pixel" at the antipodal corner. The mechanism is the symmetry structure of the basis: even-order polynomials are identical at antipodal points, odd-order polynomials are exactly inverted, and aliasing-corrupted high-order coefficients destroy the balance required to cancel intensity at the antipode. The artifact persists and grows with reconstruction order at $O$4 and $O$5 resolutions for the pixel-based method, while it is entirely absent for the edge-based method.

Corrupted QR codes. For a $O$6 QR code encoding the word "Zernike", reconstructions from pixel-based moments at $O$7 degrade to the point that the image is no longer recognized as a QR code, whereas edge-based reconstructions continue to improve with order. This is arguably the most practically consequential result in the paper: aliasing-induced reconstruction failure is not merely cosmetic but destroys machine-readable content.

Classification stability. Using 250 samples of five stroke-heavy Chinese characters under random rotation, scaling, and shearing, Zernike invariants $O$8 were evaluated via MDS embedding, ROC/AUC analysis, and nearest-centroid classification over 10,000 random hold-out partitions. Both methods perform identically up to $O$9 (accuracy ≈ 96–97%, AUC ≈ 0.83). Beyond that, the two diverge sharply:

Method	Accuracy at $e = (A, B)$0	Accuracy at $e = (A, B)$1	AUC at $e = (A, B)$2	AUC at $e = (A, B)$3
Pixel-based	96.9%	90.1%	0.83	0.78
Edge-based	96.0%	97.3%	0.83	0.84

The pixel-based invariants lose roughly 7 percentage points of accuracy as order increases, while the edge-based invariants remain stable and slightly improve. Since higher-order invariants cannot reduce the available shape information, the authors attribute the pixel-based degradation to accumulated numerical error from spatial aliasing — an interpretation consistent with the MDS projections, where class separation deteriorates only for the pixel-based variant.

Limitations and open questions

The paper is candid about several constraints. The edge-based method is one to two orders of magnitude slower than the point-mass method, so its use is justified primarily when high-order moments or exactness are required; for moderate-order analysis of high-resolution natural images, the pixel-based method remains adequate and faster. The favorable near-linear runtime depends on the empirical behavior of the adaptive quadrature (average 17 nodes), which the worst-case $e = (A, B)$4 in $e = (A, B)$5 complexity does not guarantee. The exactness claim is relative to the piecewise-constant pixel model: the method computes the exact moments of that model, not of an underlying continuous image, so any fidelity gap between the model and the true scene remains. Experiments are limited to a small number of test images and a single five-class classification task; the generality of the classification-stability result across other datasets and descriptor families is not established. Finally, the extension to 3D Zernike moments from surface meshes, and to other orthogonal moment families (Legendre, Chebyshev, Fourier–Mellin), is proposed but not carried out.

Conclusion

This paper reformulates 2D Zernike moment computation as exact boundary integration, proving that the transformed edge integrands are polynomials and therefore integrable exactly by Clenshaw–Curtis quadrature. The formulation unifies contour data, binary, grayscale, and color (quaternion) images within a single framework, and its experimental evaluation demonstrates that the two approaches coincide at low and moderate orders while the edge-based method remains stable at orders where pixel-based computation suffers aliasing-driven degradation — a 97.3% versus 90.1% classification accuracy gap at order 500, and machine-readable versus unreadable QR reconstructions. The principal trade-off is computational cost, and the principal open question is whether the boundary-integration framework extends with comparable exactness and efficiency to 3D Zernike moments and to other orthogonal moment bases.