Harmonic-based Morphometrics#

Harmonic-based morphometrics describes shape using mathematical functions that decompose outlines or surfaces into frequency components. Unlike landmark methods, harmonic approaches do not require explicit point-to-point correspondence between specimens, though specimens should represent homologous structures.

Elliptic Fourier Analysis (EFA)#

Elliptic Fourier Analysis represents closed 2D outlines as a linear combination of sine and cosine functions at various frequencies.

Mathematical Foundation#

A closed 2D outline can be parameterized by arc length \(t\) as:

\[ x(t) = A_0 + \sum_{n=1}^{N} \left( a_n \cos\frac{2\pi nt}{T} + b_n \sin\frac{2\pi nt}{T} \right) \]
\[ y(t) = C_0 + \sum_{n=1}^{N} \left( c_n \cos\frac{2\pi nt}{T} + d_n \sin\frac{2\pi nt}{T} \right) \]

where:

  • \(T\) is the total perimeter

  • \(N\) is the number of harmonics

  • \(a_n, b_n, c_n, d_n\) are Fourier coefficients

Each harmonic \(n\) contributes four coefficients that together define an ellipse. The first harmonic (\(n=1\)) captures the overall elliptical shape, while higher harmonics capture finer details.

Normalization#

Raw EFA coefficients contain information about size, rotation, and starting point. For shape analysis, coefficients are typically normalized to remove these effects:

  • Size normalization: Scale so the first harmonic has unit semi-major axis

  • Rotation normalization: Rotate so the first harmonic is aligned with a reference axis

  • Starting point normalization: Adjust phase to a standard starting point

In ktch, normalization is handled automatically:

from ktch.harmonic import EllipticFourierAnalysis

efa = EllipticFourierAnalysis(n_harmonics=20, norm=True)
coefficients = efa.fit_transform(outlines)

Choosing the Number of Harmonics#

The number of harmonics determines the level of detail captured:

Harmonics

Detail Level

Use Case

1-5

Coarse

Overall shape, major features

10-20

Moderate

Most biological applications

30+

Fine

Complex outlines, high resolution

Practical guidelines#

  • Start with 20 harmonics for most biological shapes

  • Examine reconstructed outlines to assess fit

  • More harmonics = more coefficients = higher dimensionality

  • Diminishing returns for very high harmonics

Cumulative Fourier Power#

The cumulative Fourier power indicates how much shape variance is captured by the first \(n\) harmonics. When the cumulative power exceeds 0.99, the harmonics capture 99% of the shape information.

EFA for 3D Curves#

ktch extends EFA to 3D closed curves by adding a third coordinate function \(z(t)\). This yields six coefficients per harmonic instead of four.

# 3D outline analysis
efa_3d = EllipticFourierAnalysis(n_harmonics=20, n_dim=3)
coefficients_3d = efa_3d.fit_transform(outlines_3d)

3D normalization follows the algorithm described in Godefroy et al. (2012), removing size, orientation (via Euler ZXZ decomposition), and starting point variation from the first-harmonic ellipse.

General Codomain Dimension#

The same expansion applies to data of any codomain dimension \(D\) with norm=False. Here t plays the role of the parameterization, analogous to r_theta in DHA and theta_phi in SPHARM. For \(D \in \{2, 3\}\) (spatial shape coordinates) t defaults to automatic arc-length parameterization. For \(D=1\) or \(D>3\) the codomain holds non-shape data, t is required and should be derived from the corresponding shape, onto which the data is then mapped.

Spherical Harmonic Analysis#

For 3D closed surfaces, spherical harmonics provide an analogous decomposition. A closed surface can be represented using spherical harmonic basis functions \(Y_l^m\), where \(l\) is the degree (analogous to harmonic number) and \(m\) is the order.

Usage in ktch#

ktch provides SphericalHarmonicAnalysis for working with spherical harmonic coefficients. Note that coefficient estimation requires pre-estimated surface parameterization; direct estimation from surface coordinates alone is not currently supported.

from ktch.harmonic import SphericalHarmonicAnalysis

sha = SphericalHarmonicAnalysis(n_harmonics=15)
coefficients = sha.fit_transform(parameterized_surfaces)

The codomain dimension is set by n_dim (default 3): use n_dim=1 for a scalar field on the sphere, or any \(D\) for an \(\mathbb{R}^D\)-valued field.

Registration#

SPHARM coefficients include position, orientation, and size, which will be removed for shape comparison. SphericalHarmonicAnalysis removes them through its registration parameter (e.g., "first_order" uses the degree-1 ellipsoid to canonicalize orientation, the parameter sphere, and scale). The same registration is available as a standalone transformer, SphericalHarmonicRegistration, for coefficients computed elsewhere (for example, SPHARM-PDM .coef files read with read_spharmpdm_coef):

from ktch.harmonic import SphericalHarmonicRegistration

reg = SphericalHarmonicRegistration(method="first_order", scale=False)
registered = reg.fit_transform(coefficients)

It maps coefficients to coefficients, so it composes in a scikit-learn Pipeline before PCA on a morphospace. For near-symmetric shapes, the axis and sign tie-break is ill-conditioned, so different implementations may pick different canonical frames of the same shape.

Applications#

Spherical harmonic analysis is used for:

  • Fruit shape analysis

  • Grain morphology

  • Organ shape quantification

Disk Harmonic Analysis#

For surfaces or scalar fields parameterized on a unit disk, disk harmonics provide a decomposition using Fourier–Bessel basis functions. DiskHarmonicAnalysis expands an \(\mathbb{R}^D\)-valued function on the disk, with n_dim setting the codomain dimension (n_dim=1 for a scalar field, 2/3 for planar/surface mappings). It requires a pre-estimated disk parameterization \((r, \theta)\).

from ktch.harmonic import DiskHarmonicAnalysis

dha = DiskHarmonicAnalysis(n_harmonics=10, n_dim=3)
coefficients = dha.fit_transform(surfaces, r_theta=r_theta)

Statistical Analysis#

After obtaining harmonic coefficients, standard multivariate methods apply:

Principal Component Analysis#

from sklearn.decomposition import PCA

pca = PCA(n_components=10)
pc_scores = pca.fit_transform(coefficients)

Shape Reconstruction#

Coefficients can be transformed back to outlines/surfaces for visualization:

# Reconstruct outline from coefficients
reconstructed = efa.inverse_transform(coefficients)

# Reconstruct from modified PC scores
modified_coef = pca.inverse_transform(modified_scores)
reconstructed_shape = efa.inverse_transform(modified_coef)

ktch provides built-in plot functions for these visualizations; see Morphometric Visualization for details.

Limitations#

  • Global description may miss local features

  • Sensitive to outline/surface quality and sampling

  • Normalization choices affect results

See also

References#

  • Kuhl, F. P., & Giardina, C. R. (1982). Elliptic Fourier features of a closed contour. Computer Graphics and Image Processing, 18(3), 236-258.

  • Crampton, J. S. (1995). Elliptic Fourier shape analysis of fossil bivalves. Lethaia, 28(2), 147-158.

  • Godefroy, J. E., Bornert, F., Gros, C. I., & Constantinesco, A. (2012). Elliptical Fourier descriptors for contours in three dimensions: A new tool for morphometrical analysis in biology. Comptes Rendus Biologies, 335(3), 205-213. https://doi.org/10.1016/j.crvi.2011.12.004

  • Shen, L., & Makedon, F. (2006). Spherical mapping for processing of 3D closed surfaces. Image and Vision Computing, 24(7), 743-761.

  • Verrall, S. C., & Kakarala, R. (1998). Disk-harmonic coefficients for invariant pattern recognition. Journal of the Optical Society of America A, 15(2), 389.

  • Shaqfa, M., Choi, G. P. T., Anciaux, G., & Beyer, K. (2025). Disk harmonics for analysing curved and flat self-affine rough surfaces and the topological reconstruction of open surfaces. Journal of Computational Physics, 522, 113578.