Lesson 5: Image Moments¶
Once we can isolate a blob (Lesson 4), moments let us summarize its shape with a handful of numbers: its area, centroid, orientation, and even a description that stays the same under translation, scale, and rotation. Moments offer a classic, lightweight alternative to learned features (covered later) for simple shape matching.
import numpy as np
import cv2
import matplotlib.pyplot as plt
What is a moment?¶
For a binary image, the raw moment $m_{pq}$ is defined as
$$m_{pq} = \sum_{x,y} x^p y^q \, I(x, y)$$
where $I(x,y)$ is 1 inside the shape and 0 elsewhere. A few special cases are already familiar quantities:
- $m_{00}$ = area (pixel count)
- ($\bar{x}, \bar{y}) = (m_{10}/m_{00}, m_{01}/m_{00})$ yields the centroid
The order of a moment is $p+q$. The zeroth-order moment is $m_{00}$; the first-order moments are $m_{10}$ and $m_{01}$; and the second-order moments are $m_{20}$, $m_{02}$, and $m_{11}$. cv2.moments computes all of these in one call, along with the central moments $\mu_{ij}$, which are translation-invariant; and the Hu moments, which are translation-, scale-, and rotation-invariant.
Raw moments, computed by hand¶
The formula above is just a weighted sum over pixel coordinates, so for a binary image it's nothing more than: find every foreground pixel's $(x,y)$ coordinates, then sum $x^py^q$ over them. Here it is on a tiny 3x3 image with 3 foreground pixels, computed directly with NumPy and cross-checked against cv2.moments.
im_tiny = np.array([
[0, 1, 0],
[1, 1, 0],
[0, 0, 0],
], dtype=np.uint8)
ys, xs = np.nonzero(im_tiny) # (row, col) = (y, x) of every foreground pixel
m00 = len(xs) # area: just a count of foreground pixels
m10 = xs.sum() # sum of x^1 y^0 over foreground pixels
m01 = ys.sum() # sum of x^0 y^1 over foreground pixels
m = cv2.moments(im_tiny, binaryImage=True)
print(f'foreground pixels (x, y): {[(int(x), int(y)) for x, y in zip(xs, ys)]}')
print(f'moments: m00 = {m00} m10 = {m10} m01 = {m01}')
print(f"cv2.moments: m00 = {m['m00']:.0f} m10 = {m['m10']:.0f} m01 = {m['m01']:.0f} <-- same as previous line")
print(f'centroid = ({m10 / m00:.2f}, {m01 / m00:.2f})')
An elongated, rotated blob¶
Let's draw a rotated ellipse so its orientation is easy to eyeball and check against what the moments compute.
im_binary = np.zeros((200, 200), dtype=np.uint8)
pixval = 1 # or 255, but note that it will make the moments bigger
center = (100, 100)
axes_len = (70, 25)
angle_deg = 30
cv2.ellipse(im_binary, center, axes_len, angle_deg, 0, 360, pixval, cv2.FILLED)
plt.imshow(255 * im_binary, cmap='gray')
plt.title(f'Ellipse drawn at {angle_deg} degrees')
plt.axis('off')
plt.show()
OpenCV's cv2.ellipse function produces strange effects around the border. If you want to fix these, it is easy enough to write the code to draw an ellipse yourself.
def draw_ellipse(im, cen, axes, angle_deg, val):
h, w = im.shape
Y, X = np.mgrid[0:h, 0:w]
ang = np.deg2rad(angle_deg)
XX = np.cos(ang) * (X-cen[0]) + np.sin(ang) * (Y-cen[1])
YY = -np.sin(ang) * (X-cen[0]) + np.cos(ang) * (Y-cen[1])
return np.uint8(val * (np.sqrt( (XX**2) + ((axes[0] / axes[1])**2)*(YY**2)) < (1.006)*axes[0]))
im_binary2 = np.zeros((200, 200), dtype=np.uint8)
im_binary2 = draw_ellipse(im_binary2, center, axes_len, angle_deg, pixval)
im_binary = im_binary2 # Let's replace the ellipse with our fixed ellipse
plt.imshow(255 * im_binary2, cmap='gray')
ccen = (165,90) # location of red circle when ellipse angle is 0 degrees
ccen_rotx = center[0] + np.cos(np.deg2rad(angle_deg)) * (ccen[0]-center[0]) - np.sin(np.deg2rad(angle_deg)) * (ccen[1]-center[1])
ccen_roty = center[1] + np.sin(np.deg2rad(angle_deg)) * (ccen[0]-center[0]) + np.cos(np.deg2rad(angle_deg)) * (ccen[1]-center[1])
circle = plt.Circle((ccen_rotx, ccen_roty), 10, color='red', fill=False, linewidth=1, alpha=1.0)
plt.gca().add_patch(circle)
plt.title(f'Fixed ellipse - Look closely inside red circle')
plt.axis('off')
plt.show()
Computing raw moments¶
Computing the raw moments is easy: just sum over the image pixels (as explained above).
def compute_moments(im):
summ = 0
sumx = 0
sumy = 0
for y in range(im.shape[0]):
for x in range(im.shape[1]):
v = im[y,x]
if v > 0:
summ += 1
sumx += x
sumy += y
out = {}
out['m00'] = summ
out['m10'] = sumx
out['m01'] = sumy
return out
mm = compute_moments(im_binary)
m = cv2.moments(im_binary, binaryImage=True)
print(f"moments: m00 = {mm['m00']:.0f} m10 = {mm['m10']:.0f} m01 = {mm['m01']:.0f}")
print(f"cv2.moments: m00 = {m['m00']:.0f} m10 = {m['m10']:.0f} m01 = {m['m01']:.0f} <-- same as previous line")
Area and centroid from moments¶
m = cv2.moments(im_binary, binaryImage=True)
area = m['m00']
cx = m['m10'] / m['m00']
cy = m['m01'] / m['m00']
print(f'area (m00) = {area:.0f} pixels')
print(f'centroid = ({cx:.1f}, {cy:.1f})')
plt.imshow(im_binary, cmap='gray')
plt.scatter(cx, cy, c='red', marker='x', s=80)
plt.title('Centroid from moments')
plt.axis('off')
plt.show()
Central moments¶
Raw moments $m_{pq}$ depend on where the shape happens to sit in the image — move the shape and every $m_{ij}$ (except $m_{00}$) changes. Central moments fix this by measuring around the shape's own centroid $(\bar x, \bar y)$ instead of the image origin:
$$\mu_{pq} = \sum_{x,y} (x-\bar x)^p\,(y-\bar y)^q\, I(x,y)$$
Recomputing this sum from scratch would mean re-visiting every pixel again. Instead, there's a shift-of-origin identity — the image-moment analogue of $\mathrm{Var}(X) = E[X^2] - E[X]^2$ — that computes the central moments directly from the raw moments:
$$\mu_{20} = M_{20} - \bar x\,m_{10}, \qquad \mu_{02} = m_{02} - \bar y\,m_{01}, \qquad \mu_{11} = m_{11} - \bar x\,m_{01}$$
cv2.moments returns both raw (m['m20'], ...) and central (m['mu20'], ...) moments in the same dictionary, so in practice you never need to compute these by hand —.
mu20_manual = m['m20'] - cx * m['m10']
mu02_manual = m['m02'] - cy * m['m01']
mu11_manual = m['m11'] - cx * m['m01']
print(f'moments: mu20 = {mu20_manual:.1f} mu02 = {mu02_manual:.1f} mu11 = {mu11_manual:.1f}')
print(f"cv2.moments: mu20 = {m['mu20']:.1f} mu02 = {m['mu02']:.1f} mu11 = {m['mu11']:.1f}")
Orientation from central moments¶
The central moments $\mu_{20}$, $\mu_{02}$, $\mu_{11}$ describe the spread of the shape around its centroid — essentially its covariance matrix (which we revisit in Lesson 6). The angle of the major axis (the direction of greatest spread) is
$$\theta = \frac{1}{2}\,\mathrm{atan2}\!\left(2\mu_{11},\; \mu_{20} - \mu_{02}\right)$$
theta = 0.5 * np.arctan2(2 * m['mu11'], m['mu20'] - m['mu02'])
theta_deg = np.degrees(theta)
print(f'orientation from moments = {theta_deg:.1f} degrees')
print(f'angle used to draw the ellipse = {angle_deg} degrees')
# The eigenvalues of the (normalized) covariance matrix of central moments
# give the axis lengths of the equivalent ellipse: axis = 4*sqrt(eigenvalue).
cov = np.array([[m['mu20'], m['mu11']], [m['mu11'], m['mu02']]]) / m['m00']
eigvals, _ = np.linalg.eigh(cov)
semi_major = 2 * np.sqrt(eigvals[-1])
theta = np.deg2rad(angle_deg) ####################
dx, dy = semi_major * np.cos(theta), semi_major * np.sin(theta)
print(semi_major, cx, cy, dx, dy) ######################
plt.imshow(im_binary, cmap='gray')
plt.plot([cx - dx, cx + dx], [cy - dy, cy + dy], c='red', linewidth=2)
plt.scatter(cx, cy, c='red', marker='x', s=80)
plt.title('Major axis recovered from moments')
plt.axis('off')
plt.axis('equal')
plt.show()
Hu moments: a shape descriptor invariant to pose¶
Hu moments (cv2.HuMoments) are 7 values calculated from the central moments that stay (nearly) the same regardless of the shape's position, size, and rotation. This makes them useful for comparing two shapes without first aligning them.
To see this, we build three versions of a dog shape: at the original pose, translated, and rotated + scaled. We compare the Hu moments of these versions of the shape with each other, and with a different (cat) shape.
def warp(im, mat):
rows, cols = im.shape[:2]
return cv2.warpAffine(im, np.float32(mat), (cols, rows))
im_original = cv2.imread('../img/dog_clipart.png', cv2.IMREAD_GRAYSCALE)
im_translated = warp(im_original, [[1, 0, 10], [0, 1, 20]])
im_rotated_scaled = warp(im_original, cv2.getRotationMatrix2D([60, 60], angle=45, scale=0.5))
im_different_shape = cv2.imread('../img/cat_clipart.png', cv2.IMREAD_GRAYSCALE)
fig, axes = plt.subplots(1, 4, figsize=(12, 3))
for ax, img, title in zip(
axes,
[im_original, im_translated, im_rotated_scaled, im_different_shape],
['Original', 'Translated', 'Rotated + scaled', 'Different shape'],
):
ax.imshow(img, cmap='gray')
ax.set_title(title, fontsize=10)
ax.axis('off')
plt.tight_layout()
plt.show()
Image source: Public domain pictures
Now let's compute the Hu moments.
def hu_log(im):
m = cv2.moments(im, binaryImage=True)
hu = cv2.HuMoments(m).flatten()
# log-scale since raw Hu moments span many orders of magnitude
return -np.sign(hu) * np.log10(np.abs(hu) + 1e-30)
# Format with a fixed width (":8.3f") to align the values across rows.
label_width, col_width, n_hu = 18, 8, 7
print(' ' * (label_width + 2) + 'Hu moments'.center(col_width * n_hu))
print(f'{"image":>{label_width}} ' + ''.join(f'{i:>{col_width}}' for i in range(1, n_hu + 1)))
for im, name in [
(im_original, 'original'),
(im_translated, 'translated'),
(im_rotated_scaled, 'rotated+scaled'),
(im_different_shape, 'different shape'),
]:
row = ''.join(f'{v:{col_width}.3f}' for v in hu_log(im))
print(f'{name:>{label_width}} {row}')
Note that translated and rotated+scaled versions yield similar Hu moments, whereas the different shape yields much different values. This is exactly the invariance property that makes Hu moments useful for shape matching.
Exercise¶
- Use
cv2.findContoursto get the outline of a blob from Lesson 4's binary image, then callcv2.momentson the contour instead of the full binary mask. Compare the centroid to the one computed here. - Draw a shape that is mirror-flipped rather than rotated. Are its Hu moments still close to the original? (Hint: think about what determinant/parity information moments do or do not capture.)