Lesson 6: Eigenvectors and Eigenvalues of the Moment Covariance Matrix¶
In Lesson 5 we used central moments to find a blob's orientation and, with a bit of extra algebra, its semi-axis lengths. This lesson makes that connection explicit and general: the central moments of any blob (not just an ellipse) define a $2\times2$ covariance matrix, and the eigenvectors/eigenvalues of that matrix directly give the orientation and size of the equivalent ellipse — the ellipse with the same area and same second-moment spread as the blob.
Refresher: what are eigenvectors and eigenvalues?¶
For a square matrix $A$, a nonzero vector $v$ is an eigenvector with eigenvalue $\lambda$ if
$$Av = \lambda v$$
In words: applying $A$ to $v$ doesn't rotate $v$ off its own line — it only scales it, by a factor of $\lambda$. Most vectors get both rotated and scaled by $A$; eigenvectors are the special directions that only get scaled.
For a $2\times2$ symmetric matrix (like the covariance matrices we'll build below), something even nicer is guaranteed: there are always two eigenvectors, they are perpendicular to each other, and their eigenvalues are real numbers. Geometrically, $A$ takes a circle of unit vectors and stretches it into an ellipse whose axes point along the eigenvectors, with each semi-axis length equal to the corresponding eigenvalue. That picture — eigenvectors as axis directions, eigenvalues as axis lengths — is exactly what we'll exploit for shape analysis below.
import numpy as np
import cv2
import matplotlib.pyplot as plt
A = np.array([[3.0, 1.0],
[1.0, 1.5]])
eigvals, eigvecs = np.linalg.eigh(A) # ascending order
print('eigenvalues :', eigvals)
print('eigenvectors (columns):\n', eigvecs)
theta = np.linspace(0, 2 * np.pi, 200)
circle = np.stack([np.cos(theta), np.sin(theta)]) # unit circle, as column vectors
ellipse = A @ circle # apply A to every point on the circle
fig, ax = plt.subplots(figsize=(4.5, 4.5))
ax.plot(circle[0], circle[1], '--', color='gray', label='unit circle')
ax.plot(ellipse[0], ellipse[1], color='#3498db', label='A @ circle')
for val, vec in zip(eigvals, eigvecs.T):
ax.plot([0, val * vec[0]], [0, val * vec[1]], color='red', linewidth=2)
ax.scatter(0, 0, color='black', zorder=5)
ax.set_aspect('equal')
ax.legend(loc='upper left', fontsize=8)
ax.set_title('A circle stretched by A; red lines = eigenvectors x eigenvalues')
plt.show()
From one eigenvector to both at once: diagonalization¶
$Av = \lambda v$ above holds for a single eigenvector. Stack both eigenvectors as the columns of a matrix $P = \begin{bmatrix}v_1 & v_2\end{bmatrix}$, and both eigenvalues on the diagonal of $\Lambda = \begin{bmatrix}\lambda_1 & 0\\0 & \lambda_2\end{bmatrix}$. Then $Av_1=\lambda_1 v_1$ and $Av_2=\lambda_2 v_2$, side by side, become one matrix equation:
$$AP = P\Lambda$$
Because $A$ is symmetric, its eigenvectors are not just perpendicular but orthonormal, which makes $P$ an orthogonal matrix (inverting it is just a transpose: $P^{-1} = P^\top$). That lets us solve for $A$ itself:
$$A = P\Lambda P^\top$$
This is the diagonalization of $A$, and geometrically it's three steps: $P^\top$ rotates coordinates into the eigenvector frame, $\Lambda$ scales along those (now axis-aligned) directions independently, and $P$ rotates back. Read in that rotated frame, $A$ is genuinely diagonal — the off-diagonal entry, which mixes $x$ and $y$ together, is exactly zero. When $A$ is a covariance matrix, that off-diagonal entry is $\mu_{11}$: diagonalizing the covariance matrix and finding the ellipse's own natural axes are the same operation — both mean rotating into the one coordinate frame where the blob's spread in $x$ and $y$ no longer mixes at all.
P = eigvecs
Lam = np.diag(eigvals)
print('A @ P:\n', np.round(A @ P, 4))
print('P @ Lambda:\n', np.round(P @ Lam, 4))
print('max |A@P - P@Lambda| (should be ~0):', np.abs(A @ P - P @ Lam).max())
print()
reconstructed = P @ Lam @ P.T
print('P @ Lambda @ P.T (reconstructed A):\n', np.round(reconstructed, 4))
print('original A:\n', A)
print('max |reconstructed - A|:', np.abs(reconstructed - A).max())
print()
diagonalized = P.T @ A @ P # A, read in the eigenvector frame
print('P.T @ A @ P (A in its own eigenbasis -- should be diagonal):\n', np.round(diagonalized, 4))
From central moments to a covariance matrix¶
Treat a blob's pixels as samples from a 2D distribution. Its covariance matrix, in terms of the central moments $\mu_{ij}$ and area $\mu_{00}=m_{00}$, is
$$A = \frac{1}{m_{00}}\begin{bmatrix}\mu_{20} & \mu_{11} \\ \mu_{11} & \mu_{02}\end{bmatrix}$$
This is exactly the same formula used for the covariance matrix of a scatter of points $(x,y)$, just weighted by pixel membership instead of by sample index. The off-diagonal $\mu_{11}$ term is precisely the entry that mixing $x$ and $y$ leaves behind, so this is exactly the matrix we want to diagonalize: its eigenvectors will point along the blob's own major and minor axes — the directions of greatest and least spread — and its eigenvalues will measure how much spread there is along each.
def covariance_from_moments(binary):
m = cv2.moments(binary, binaryImage=True)
cx, cy = m['m10'] / m['m00'], m['m01'] / m['m00']
cov = np.array([[m['mu20'], m['mu11']],
[m['mu11'], m['mu02']]]) / m['m00']
return cov, (cx, cy)
Eigen decomposition gives orientation and size¶
np.linalg.eigh returns eigenvalues in ascending order along with their eigenvectors (as columns). For a filled ellipse with semi-axes $a \ge b$, the eigenvalues work out to $\lambda_{\max} = a^2/4$ and $\lambda_{\min}=b^2/4$, so
$$a = 2\sqrt{\lambda_{\max}}, \qquad b = 2\sqrt{\lambda_{\min}}$$
and the corresponding eigenvectors point along the major and minor axes. Note that the angle of the major axis is the same as the angle computed in Lesson 5.
def principal_axes(binary):
cov, center = covariance_from_moments(binary)
eigvals, eigvecs = np.linalg.eigh(cov) # ascending order
semi_axes = 2 * np.sqrt(np.clip(eigvals, 0, None))
# reorder so index 0 is major (largest), index 1 is minor
order = [1, 0]
# eigvecs[:, order] has the two eigenvectors as its COLUMNS; transpose so that
# unpacking it below (which iterates over ROWS) hands back the two actual eigenvectors,
# not a scrambled mix of their x- and y-components
return center, semi_axes[order], eigvecs[:, order].T
def draw_axes(ax, binary, color='red'):
center, (a, b), (v_major, v_minor) = principal_axes(binary)
cx, cy = center
ax.imshow(binary, cmap='gray')
for length, vec, lw in [(a, v_major, 2.5), (b, v_minor, 1.5)]:
dx, dy = length * vec
ax.plot([cx - dx, cx + dx], [cy - dy, cy + dy], c=color, linewidth=lw)
ax.scatter(cx, cy, c=color, marker='x', s=60)
ax.axis('off')
return center, (a, b)
Sanity check: a known ellipse¶
As in Lesson 5, we draw an ellipse with known parameters and confirm the eigen-based estimate recovers them — but this time we plot both axes, not just the major one.
binary = np.zeros((200, 200), dtype=np.uint8)
cv2.ellipse(binary, (100, 100), (70, 25), 30, 0, 360, 255, -1)
fig, ax = plt.subplots(figsize=(4, 4))
center, (a, b) = draw_axes(ax, binary)
ax.set_title('Major (thick) and minor (thin) axes from eigenvectors')
plt.show()
print(f'recovered semi-axes: a={a:.1f}, b={b:.1f} (drawn with 70, 25)')
It works on arbitrary shapes, too¶
The real power of this approach is that it doesn't require the blob to be an ellipse at all. Every binary shape has some equivalent ellipse — the one that matches its area, centroid, and second-moment spread. This gives a compact 5-number summary (center, two semi-axis lengths, orientation) of an arbitrarily shaped blob.
im_glasses = cv2.imread('../img/glasses_outline.png', cv2.IMREAD_GRAYSCALE)
shapes = {'Glasses': im_glasses}
fig, ax = plt.subplots(figsize=(5, 4))
center, (a, b) = draw_axes(ax, im_glasses)
ax.set_title(f'Glasses with axes\na={a:.0f}, b={b:.0f}', fontsize=10)
plt.show()
Notice the equivalent ellipse doesn't try to trace the shape's boundary — it summarizes the distribution of mass around the centroid, so it cuts across the narrow bridge between the two lenses rather than following it, and its outline falls outside the lenses near the top and bottom where the glasses taper.
Eccentricity: how elongated is a shape?¶
The ratio of the eigenvalues (or semi-axes) gives a single number describing elongation, independent of orientation and overall size:
$$\text{eccentricity} = \sqrt{1 - \frac{\lambda_{\min}}{\lambda_{\max}}}$$
This ranges from 0 (a circle, both axes equal) to nearly 1 (a very thin, elongated shape).
print(f'{"shape":>10} {"a":>6} {"b":>6} {"eccentricity":>13}')
for name, img in shapes.items():
_, (a, b), _ = principal_axes(img)
ecc = np.sqrt(1 - (b / a) ** 2)
print(f'{name:>10} {a:6.1f} {b:6.1f} {ecc:13.3f}')
Exercise¶
- Draw a perfect circle and confirm its eccentricity is (close to) 0.
- Overlay the actual fitted ellipse outline (not just the axis lines) using
cv2.ellipsewith the center,(2a, 2b)as the full axes lengths, and the orientation angle from the eigenvectors. Compare it tocv2.fitEllipseapplied to the shape's contour — do they agree? - Run
principal_axeson the filled blob ofglasses(after flood fill)). How much doa,b, and the eccentricity change? Does that match your intuition for how moments depend on where mass sits, not just the overall silhouette?