Theory · Shape · Elliptic Fourier

A seed is a sum of ellipses.

Walk once around a seed outline, noting x and y. Both lists repeat on every lap, and whatever repeats becomes a Fourier series. Written with x and y kept separate, each term of the series draws an ellipse, and half a dozen ellipses already describe the shape.

Go deeper · mathematics and references

The outline as a curve

A seed outline is a closed curve. To compute with it, walk along the outline at constant speed from any starting point and call \(t\) the length covered so far. After one lap, \(t\) reaches the perimeter \(T\) and the walk starts over at the same place.

Definition

An outline parametrized by arc length is a pair of functions \(x(t)\) and \(y(t)\), with \(0 \le t \lt T\), such that the point \(\big(x(t), y(t)\big)\) moves along the outline at speed 1 and returns to the start at \(t = T\). Both functions are periodic: \(x(t+T) = x(t)\) and \(y(t+T) = y(t)\).

In practice the outline arrives as a polygon, the list of vertices that comes out of segmentation. Between two vertices, \(x\) and \(y\) vary linearly with \(t\). They are continuous functions with corners at the vertices, and it becomes clear further on that this is enough for the Fourier series to converge to them at every point.

Two Fourier series

Each coordinate has its own series. With \(\theta = 2\pi t/T\), the angle that measures the fraction of the lap already covered,

\[ x(t) = A_0 + \sum_{n=1}^{\infty} \big( a_n \cos n\theta + b_n \sin n\theta \big), \qquad y(t) = C_0 + \sum_{n=1}^{\infty} \big( c_n \cos n\theta + d_n \sin n\theta \big) \](1)

with the usual coefficients:

\[ a_n = \frac{2}{T}\int_0^T x(t)\cos n\theta\,dt, \qquad b_n = \frac{2}{T}\int_0^T x(t)\sin n\theta\,dt, \qquad A_0 = \frac{1}{T}\int_0^T x(t)\,dt \](2)

and \(c_n\), \(d_n\), \(C_0\) are the same, with \(x\) replaced by \(y\). The pair \((A_0, C_0)\) is the center of the outline, the mean of the points along the perimeter. The term of order \(n\) makes \(n\) turns while the point makes one: \(n = 1\) describes the coarse shape, and the higher-order terms describe the details of the edge.

Each term is an ellipse

Kuhl and Giardina's trick [1] is to look at term \(n\) of both series at once, as a 2 × 2 matrix applied to a point on the circle:

\[ \begin{pmatrix} x_n \\ y_n \end{pmatrix} = M_n \begin{pmatrix} \cos n\theta \\ \sin n\theta \end{pmatrix}, \qquad M_n = \begin{pmatrix} a_n & b_n \\ c_n & d_n \end{pmatrix} \](3)

Result 1 · the harmonic is an ellipse

As \(t\) goes once around the outline, the harmonic \(n\) traces an ellipse centered at the origin, \(n\) times. Its semi-axes are the singular values \(\sigma_1 \ge \sigma_2\) of \(M_n\), its area is \(\pi\,\lvert\det M_n\rvert\) and the direction of rotation is given by the sign of \(\det M_n\).

Sketch. The point \((\cos n\theta, \sin n\theta)\) moves on the unit circle. Write the singular value decomposition \(M_n = U\Sigma V^\top\): \(V^\top\) rotates the circle, \(\Sigma\) stretches the axes by \(\sigma_1\) and \(\sigma_2\), and \(U\) rotates again. The circle becomes an ellipse with semi-axes \(\sigma_1\) and \(\sigma_2\), and the area is multiplied by \(\lvert\det M_n\rvert = \sigma_1\sigma_2\). A negative determinant reverses the direction.

The whole curve is the sum of the ellipses. The center of the second ellipse moves along the first, the center of the third moves along the second, and the tip of the last one draws the outline. The instrument below draws this chain turning over a real seed.

Result 2 · the area comes from the coefficients

The area inside the outline, positive when the outline is traversed counterclockwise, is

\[ A = \pi \sum_{n=1}^{\infty} n\,\big(a_n d_n - b_n c_n\big) = \pi \sum_{n=1}^{\infty} n \det M_n \](4)

Sketch. By Green's theorem, \(A = \tfrac12\oint (x\,dy - y\,dx)\). Substitute the series (1) and integrate with respect to \(\theta\) from 0 to \(2\pi\). Products of terms of different orders integrate to zero, and each order \(n\) contributes \(\pi n (a_n d_n - b_n c_n)\).

For the seed in Figure 1, formula (4) with 40 harmonics gives 727.12 px², and the shoelace formula applied to the polygon vertices gives 727.14 px². The difference, 0.003%, is the tail of the series that was left out.

From polygon to coefficients

The integrals (2) seem to call for numerical integration. For a polygon they do not. Let the polygon have \(K\) sides, with \(\Delta x_p\) and \(\Delta y_p\) the displacements of side \(p\), \(\Delta t_p = \sqrt{\Delta x_p^2 + \Delta y_p^2}\) its length, \(t_p = \Delta t_1 + \dots + \Delta t_p\) and \(\theta_p = 2\pi t_p / T\). Kuhl and Giardina [1] arrived at

\[ a_n = \frac{T}{2\pi^2 n^2}\sum_{p=1}^{K}\frac{\Delta x_p}{\Delta t_p}\big(\cos n\theta_p - \cos n\theta_{p-1}\big), \qquad b_n = \frac{T}{2\pi^2 n^2}\sum_{p=1}^{K}\frac{\Delta x_p}{\Delta t_p}\big(\sin n\theta_p - \sin n\theta_{p-1}\big) \](5)

and \(c_n\), \(d_n\) likewise with \(\Delta y_p\) in place of \(\Delta x_p\). The center is the mean of each side weighted by its length,

\[ A_0 = \frac{1}{T}\sum_{p=1}^{K} \Delta t_p\,\frac{x_{p-1} + x_p}{2}, \qquad C_0 = \frac{1}{T}\sum_{p=1}^{K} \Delta t_p\,\frac{y_{p-1} + y_p}{2}, \](6)

which is the closed form of (2) for a function that is linear on each side, equivalent to the one in the original paper.

Where (5) comes from

Integrate (2) by parts. Since \(x\) is periodic, the boundary term vanishes and what remains is \(a_n = -\frac{1}{\pi n}\int_0^T x'(t)\sin n\theta\,dt\). On a side of the polygon, \(x'(t) = \Delta x_p/\Delta t_p\) is constant, and the integral of \(\sin n\theta\) from \(t_{p-1}\) to \(t_p\) equals \(-\frac{T}{2\pi n}\big(\cos n\theta_p - \cos n\theta_{p-1}\big)\). Summing over the sides gives (5).

Two practical consequences. The computation is exact for the polygon, with no resampling and no numerical integration, and costs \(K \cdot N\) sines and cosines. And the factor \(1/n^2\) shows that the coefficients of a polygon decay at least like \(1/n^2\). The sum of their magnitudes converges, so the series converges uniformly [9]: the reconstruction gets arbitrarily close to the polygon at every point, corners included.

Reconstructing with N harmonics

Truncating the series after \(N\) terms gives the curve \(z_N(t) = \big(x_N(t), y_N(t)\big)\), the sum of the center and the first \(N\) ellipses. It stores \(4N + 2\) numbers. The question is how much is lost.

Result 3 · truncating gives the best approximation

Among all curves whose coordinates are trigonometric polynomials of degree \(N\), the truncated series is the closest to the outline in mean squared error, and that error is the tail of the coefficients:

\[ e_N^2 = \frac{1}{T}\int_0^T \big\lvert z(t) - z_N(t) \big\rvert^2\,dt = \frac12 \sum_{n \gt N} \big( a_n^2 + b_n^2 + c_n^2 + d_n^2 \big) \](7)

Sketch. The functions \(1, \cos n\theta, \sin n\theta\) are orthogonal in \([0, T)\). Truncating the series amounts to projecting \(x\) and \(y\) orthogonally onto the space of trigonometric polynomials of degree \(N\), and the orthogonal projection is the nearest point. What is left over is orthogonal to what stays, and the Pythagorean theorem in this space (Parseval's identity) gives (7) [9].

The tail has infinitely many terms, but there is no need to sum it. The total energy \(\tfrac{1}{T}\int_0^T \lvert z - z_0 \rvert^2 dt\), with \(z_0 = (A_0, C_0)\), can be computed exactly from the polygon, one side at a time, and \(e_N^2\) is that energy minus half the sum of the squares of the first \(N\) harmonics. This is how the instrument computes the error. The same statement reappears on the SVD page under another name: there, truncating an orthogonal decomposition gives the best low-rank approximation.

Real dataFigure 1. A seed of Amorpha fruticosa, one of the 145 silhouettes, reconstructed with 1, 2, 3, 5, 10 and 20 harmonics. Dashed: the 128-point outline. Orange: the center plus the first N ellipses. The error is \(e_N\) from equation (7), as a percentage of the seed length (51.1 px).

With \(N = 1\) we get the ellipse that best summarizes the seed, and the error is 6.3% of the length. With \(N = 3\) the kidney-shaped indentation appears, and the error falls to 1.9%. With 10 harmonics the curve is within 0.48% of the outline, and with 20, within 0.21%.

Real dataFigure 2. Error \(e_N\) against the number of harmonics for the 145 silhouettes, on a log scale. The band runs from the first to the third quartile. The thin lines mark one pixel and half a pixel of the silhouette, which is about 51 px in length.

At the median of the 145 seeds, the error falls from 3.45% with one harmonic to 1.40% with three and 0.83% with five. Half of the seeds fall below 1% with \(N = 5\), and 90% of them with \(N = 10\). The most demanding one, a Calligonum mongolicum with an indented outline, only gets there with 17 harmonics, and the Bromus inermis on the dashed line, with 14.

Result 4 · a centrally symmetric shape has no even harmonics

If the outline is symmetric about the center, that is, \(z(t + T/2) - z_0 = -\big(z(t) - z_0\big)\), then \(a_n = b_n = c_n = d_n = 0\) for every even \(n\).

Sketch. Replace \(t\) with \(t + T/2\) in (2). The cosine gets the factor \(\cos(n\theta + n\pi) = (-1)^n \cos n\theta\) and the function changes sign around the center, hence \(a_n = -(-1)^n a_n\). For even \(n\) this forces \(a_n = 0\). The same holds for \(b_n\), \(c_n\) and \(d_n\).

A rice grain is nearly symmetric about its center. In the 40 rice silhouettes, the median size of the second harmonic is 1.35% of the length, and that of the third is 2.97%. That is why the rice error barely moves from \(N = 1\) to \(N = 2\) (from 3.57% to 3.23%) and plunges with \(N = 3\) (1.40%). It is the same step that the median in Figure 2 shows between 2 and 3.

Instrument · ellipses over a real outline

harmonics N 3 show the ellipses normalize by the 1st ellipse

–error \(e_N\), % of length

–stored numbers, \(4N+2\)

–area by formula (4) ÷ polygon area

–1st ellipse, minor ÷ major semi-axis

Size of each harmonic, \(\sqrt{(a_n^2+b_n^2+c_n^2+d_n^2)/2}\), as % of length (log scale)

The coefficients are computed here, in your browser, by formula (5), from the 128 points of each outline embedded in the page. Dashed is the outline; orange, the sum of the N ellipses; blue, the chained ellipses, each centered on the tip of the previous one. With "normalize", the figure is rotated, scaled and translated so that the first ellipse lies flat, with semi-major axis 1, and the moving point starts at the tip of the major axis. In soybean, the ratio of the semi-axes is close to 1 and the first ellipse almost becomes a circle.

Normalizing by the first ellipse

The same seed, rotated or with the outline starting at another point, gives different coefficients. To compare shapes, position, size, rotation and starting point have to be taken out of the computation. Kuhl and Giardina do this using only the first ellipse [1]. First, the angle \(\theta_1\) that takes the start of the outline to the tip of the major axis of the first ellipse:

\[ \theta_1 = \tfrac12 \operatorname{atan2}\!\big(2(a_1 b_1 + c_1 d_1),\; a_1^2 - b_1^2 + c_1^2 - d_1^2\big) \](8)

Then, with \(R(\varphi)\) the rotation of the plane by the angle \(\varphi\), each harmonic is shifted in time by \(n\theta_1\), rotated by the angle \(\psi_1\) of the major axis and divided by the semi-major axis \(E\):

\[ \begin{gathered} M_n^{*} = M_n\,R(n\theta_1), \qquad \psi_1 = \operatorname{atan2}(c_1^{*}, a_1^{*}), \qquad E = \sqrt{(a_1^{*})^2 + (c_1^{*})^2}, \\[4pt] M_n'' = \frac{1}{E}\,R(-\psi_1)\,M_n^{*} \end{gathered} \](9)

Result 5 · normalizing is a 2 × 2 SVD

After (9), the first matrix becomes \(M_1'' = \begin{pmatrix} 1 & 0 \\ 0 & \pm\sigma_2/\sigma_1 \end{pmatrix}\): the first ellipse lies flat along the \(x\) axis, with semi-major axis 1, and the outline starts at the tip of that axis.

Sketch. Write \(M_1 = U\Sigma V^\top\). The point of the first ellipse farthest from the center corresponds to the direction \(v_1\), the first column of \(V\), and (8) is the angle of \(v_1\), the leading eigenvector of \(M_1^\top M_1\). Then \(M_1 R(\theta_1)\) has first column \(\sigma_1 u_1\), \(\psi_1\) is the angle of \(u_1\) and \(E = \sigma_1\). Rotating by \(-\psi_1\) and dividing by \(\sigma_1\) leaves \(\operatorname{diag}(1, \pm\sigma_2/\sigma_1)\). It is the same decomposition as on the next page, in miniature.

What is lost is accounted for. Out go the position \((A_0, C_0)\), the size \(E\), the rotation \(\psi_1\) and the starting point \(\theta_1\), and three numbers of the first ellipse are fixed (\(a_1'' = 1\), \(b_1'' = c_1'' = 0\)). Of the \(4N + 2\) numbers, \(4N - 3\) remain, describing only the shape. Two ambiguities are left. Swapping \(\theta_1\) for \(\theta_1 + \pi\) also takes the start to an end of the major axis, and the two choices differ by the sign of the even harmonics. And a seed that lands with the other face up appears mirrored, with different coefficients: the normalization does not recognize the mirror image as the same shape.

We checked the invariance with the seed in Figure 1. After rotating the outline by 40°, enlarging it 1.7 times and starting the count a third of the way around the lap, the normalized coefficients changed by at most 6 × 10⁻¹⁶, the size of machine rounding error. The first three normalized harmonics are these:

Normalized coefficients of the *Amorpha fruticosa* in Figure 1
n	\(a_n''\)	\(b_n''\)	\(c_n''\)	\(d_n''\)
1	1	0	0	0,490
2	0,024	−0,049	−0,138	0,036
3	0,086	0,013	0,089	0,026

The \(d_1'' = 0{,}490\) is the ratio between the semi-axes of the first ellipse. The second harmonic, which would not exist in a centrally symmetric shape, has \(c_2'' = -0{,}138\): it is the asymmetry between the concave and convex sides of the kidney.

Where it fails · nearly round seeds

The normalization depends on the direction of the major axis of the first ellipse, and that direction is well defined only when \(\sigma_1\) is much larger than \(\sigma_2\). A small perturbation in the matrix can rotate the singular vectors by an angle of the order of the perturbation divided by \(\sigma_1 - \sigma_2\) [10]. In 24 of the 33 soybean silhouettes, the minor semi-axis of the first ellipse exceeds 90% of the major (median 0.928). Changing the blur used to extract the outline from 0.6 to 1.0 px rotates the first ellipse of the soybean seeds by 1.6° at the median and by up to 9.0°. In rice, the same change rotates it by 0.05°, and at most 0.2°. In soybean, the normalized coefficients carry this rotation. For nearly round seeds, compare quantities that do not depend on the angle, such as the semi-axes \(\sigma_1\) and \(\sigma_2\) of each harmonic, or normalize by a reference marked on the seed.

The classical descriptors

Before 1982 there were two families of Fourier descriptors. Zahn and Roskies [2] expanded the tangent angle along the outline as a series, a function that does not change under translation or a change of scale and that rotation merely shifts by a constant. The price is a dependence on the direction of the tangent, which is sensitive to the pixel staircase, and the truncated series does not guarantee a closed curve. Granlund [3] and Persoon and Fu [4] treated the outline as a complex number \(z = x + iy\), with a single series:

\[ z(t) = \sum_{k=-\infty}^{\infty} Z_k\, e^{ik\theta} \](10)

Each term is a circle of radius \(\lvert Z_k \rvert\) that turns \(k\) times per lap, counterclockwise for \(k \gt 0\) and clockwise for \(k \lt 0\). This is what the video shows, on a model outline.

Mathematical animationVideo. The complex series (10) on a model outline, with 3, 5, 7, 11 and 25 terms. Not a photographed seed. 57 s.

Result 6 · an ellipse is two circles

The elliptic harmonic \(n\) and the pair of complex terms \(\pm n\) carry the same information:

\[ Z_n = \tfrac12\big[(a_n + d_n) + i\,(c_n - b_n)\big], \qquad Z_{-n} = \tfrac12\big[(a_n - d_n) + i\,(c_n + b_n)\big] \](11)

and the semi-axes of the ellipse are \(\sigma_1 = \lvert Z_n\rvert + \lvert Z_{-n}\rvert\) and \(\sigma_2 = \big\lvert \lvert Z_n\rvert - \lvert Z_{-n}\rvert \big\rvert\).

Sketch. In (3), write \(\cos n\theta = (e^{in\theta} + e^{-in\theta})/2\) and \(\sin n\theta = (e^{in\theta} - e^{-in\theta})/2i\) and group the terms in \(e^{\pm in\theta}\). For the semi-axes, check that \(\lvert Z_n\rvert^2 - \lvert Z_{-n}\rvert^2 = a_n d_n - b_n c_n = \det M_n\) and that \(2(\lvert Z_n\rvert^2 + \lvert Z_{-n}\rvert^2) = a_n^2 + b_n^2 + c_n^2 + d_n^2 = \sigma_1^2 + \sigma_2^2\).

The three terms at the start of the video (\(k = -1, 0, 1\)) are the center and one ellipse. The 25 at the end are the center and 12 ellipses. The elliptic formulation separates what each axis does and gives a normalization with geometric meaning, and it is the one that spread in plant morphometry. Rohlf and Archie [5] compared Fourier methods on mosquito wings and found the elliptic descriptors the most promising. The SHAPE program [6] and the Momocs package [7] perform the computation shown on this page directly from images. In wheat seeds, elliptic descriptors found more regions of the genome linked to grain shape than the length-to-width ratio did [8].

In SeedCounter

The machine proposes the outline of each seed and the person checks it. The checked outline is a polygon, and the polygon is the exact input of formula (5): the coefficients come out of it with no resampling and no numerical integration. The app measures each seed on a simplified outline of up to 48 sides. With 48 vertices there are 48 values of \(x\) and 48 of \(y\), and each harmonic uses two of each. From order 24 on, the series starts describing the corners of the polygon, not the seed. Elliptic Fourier is not in the app yet, either as a computation or as a spreadsheet column. This page shows the computation the app would run and how many harmonics make sense at each resolution.

Where it fails · Fourier copies the outline it receives

The series faithfully describes the shape it is given, including the wrong one. If a shadow stuck to the edge, or two touching seeds became one, the coefficients describe the shadow and the pair. The post Light decides the outline shows where much of this error comes from, which is why checking comes before the computation.

Resolution matters too. In the silhouettes on this page, about 51 px long, the median error stays below half a pixel from \(N = 5\) on. What the following harmonics add is smaller than the image resolution and ends up describing the pixel staircase and the blur. And the series compares points by arc length, not by homologous points: the hilum of two seeds can fall at different values of \(t\).

Data

The 145 silhouettes of 64 × 64 px from the public SeedCounter examples (rice, soybean, coffee, corn and native species), the same ones as on the Learn page. Each silhouette is laid flat along its major axis and scaled to about 51 px long, so the measurements are in silhouette px, not mm. The outline is taken from the largest region of each silhouette, with holes filled, after a Gaussian blur of 0.6 px, and resampled to 128 equally spaced points. The length is the maximum Feret diameter. The script that regenerates the figures is in _src/figuras/teoria-forma-fourier.py.

References

Kuhl F.P., Giardina C.R. (1982). Elliptic Fourier features of a closed contour. Computer Graphics and Image Processing 18(3), 236–258. doi:10.1016/0146-664X(82)90034-X
Zahn C.T., Roskies R.Z. (1972). Fourier descriptors for plane closed curves. IEEE Transactions on Computers C-21(3), 269–281. doi:10.1109/TC.1972.5008949
Granlund G.H. (1972). Fourier preprocessing for hand print character recognition. IEEE Transactions on Computers C-21(2), 195–201. doi:10.1109/TC.1972.5008926
Persoon E., Fu K.S. (1977). Shape discrimination using Fourier descriptors. IEEE Transactions on Systems, Man, and Cybernetics 7(3), 170–179. doi:10.1109/TSMC.1977.4309681
Rohlf F.J., Archie J.W. (1984). A comparison of Fourier methods for the description of wing shape in mosquitoes (Diptera: Culicidae). Systematic Zoology 33(3), 302–317. doi:10.2307/2413076
Iwata H., Ukai Y. (2002). SHAPE: a computer program package for quantitative evaluation of biological shapes based on elliptic Fourier descriptors. Journal of Heredity 93(5), 384–385. doi:10.1093/jhered/93.5.384
Bonhomme V., Picq S., Gaucherel C., Claude J. (2014). Momocs: outline analysis using R. Journal of Statistical Software 56(13). doi:10.18637/jss.v056.i13
Williams K., Munkvold J., Sorrells M. (2013). Comparison of digital image analysis using elliptic Fourier descriptors and major dimensions to phenotype seed shape in hexaploid wheat (Triticum aestivum L.). Euphytica 190(1), 99–116. doi:10.1007/s10681-012-0783-0
Stein E.M., Shakarchi R. (2003). Fourier Analysis: An Introduction. Princeton University Press.
Golub G.H., Van Loan C.F. (2013). Matrix Computations, 4th ed. Johns Hopkins University Press.

First the outline, then the ellipses.

SeedCounter proposes the outline of each seed and you check it. Only after that does the shape become a number.

Open SeedCounter