Lesson 26: Epipolar Geometry¶
Lesson 21's stereo matching assumed a rectified pair, where corresponding points always fall on the same row — a special case. This lesson covers the general two-view relationship between any pair of images of a static scene, calibrated or not: the fundamental matrix and essential matrix, which say that a point's match in the other image isn't just somewhere — it's constrained to lie on a particular line.
import numpy as np
import cv2
import matplotlib.pyplot as plt
The epipolar constraint¶
For a point $x_1$ in image 1 and its true match $x_2$ in image 2 (both in homogeneous pixel coordinates), there's a $3\times3$ matrix $F$ — the fundamental matrix — such that
$$x_2^\top F x_1 = 0$$
for every corresponding pair, regardless of scene geometry. $F$ depends only on the two cameras' relative pose and (uncalibrated) intrinsics. Rearranged, $l_2 = F x_1$ is the epipolar line in image 2: the 1D line along which $x_1$'s match is guaranteed to lie. This is exactly the mechanism that made Lesson 21's row-restricted search valid — rectification is nothing more than choosing $F$ to be a pure horizontal-line generator.
A synthetic two-camera scene¶
We build two cameras with known intrinsics $K$ and a known relative rotation/translation, project the same random 3D points into both, and use the resulting correspondences to recover $F$ — so we always have ground truth to check against.
rng = np.random.default_rng(0)
K = np.array([[500, 0, 320], [0, 500, 240], [0, 0, 1]], dtype=np.float64)
R1, t1 = np.eye(3), np.zeros(3) # camera 1: at the origin, looking down +z
angle = np.radians(15)
R_true = np.array([[np.cos(angle), 0, np.sin(angle)],
[0, 1, 0],
[-np.sin(angle), 0, np.cos(angle)]])
t_true = np.array([0.5, 0.0, 0.1]) # camera 2: rotated 15 deg, shifted along x
P1 = K @ np.hstack([R1, t1.reshape(3, 1)])
P2 = K @ np.hstack([R_true, t_true.reshape(3, 1)])
def project(P, points_3d):
homogeneous = np.hstack([points_3d, np.ones((len(points_3d), 1))])
projected = (P @ homogeneous.T).T
return projected[:, :2] / projected[:, 2:3]
points_3d = rng.uniform(-1, 1, (40, 3)) + np.array([0, 0, 5]) # in front of both cameras
x1 = project(P1, points_3d)
x2 = project(P2, points_3d)
The 8-point algorithm¶
Each correspondence gives one linear equation in the 9 unknown entries of $F$ (expanding $x_2^\top F x_1 = 0$). With 8 or more correspondences, this becomes a homogeneous least-squares problem solvable by SVD — the same kind of DLT machinery used for the homography fit in Lesson 23, just with a different constraint equation and an extra step to enforce that $F$ has rank 2 (since $F$ is singular by construction — a fundamental matrix always has a zero eigenvalue).
The SVD, briefly¶
Any matrix $A$ factors as $A=U\Sigma V^\top$, with $U$ and $V$ orthogonal and $\Sigma$ diagonal, non-negative, sorted largest to smallest (the singular values). Two facts about this are all the code below actually needs:
- Solving $Ax=0$ as closely as possible, subject to $\|x\|=1$: the answer is $V$'s last column (equivalently $V^\top$'s last row) — the right singular vector with the smallest singular value. This is Lesson 22's total-least-squares trick, generalized: $V$'s columns are the eigenvectors of $A^\top A$, sorted by eigenvalue, so the smallest-singular-value direction is exactly the same "direction of least fit error" as the smallest-eigenvalue eigenvector of a small covariance matrix — just for a $9\times9$ one here, which nobody wants to form and eigendecompose by hand.
- Enforcing a rank constraint: zeroing the smallest singular value(s) and reconstructing gives the closest lower-rank matrix (in a least-squares sense) to the original. That's exactly what's needed to force the estimated $F$ — which must be exactly rank 2 — back onto that constraint after the unconstrained 9-parameter solve inevitably drifts off it slightly.
def normalize_points(x):
"""Shift/scale so points are centered at the origin with average distance sqrt(2) -- standard
numerical-conditioning trick (Hartley normalization) for the 8-point algorithm."""
mean = x.mean(axis=0)
std = x.std()
T = np.array([[1 / std, 0, -mean[0] / std],
[0, 1 / std, -mean[1] / std],
[0, 0, 1]])
x_h = np.hstack([x, np.ones((len(x), 1))])
return (T @ x_h.T).T, T
def eight_point_algorithm(x1, x2):
x1n, T1 = normalize_points(x1)
x2n, T2 = normalize_points(x2)
A = np.array([[xb * xa, xb * ya, xb, yb * xa, yb * ya, yb, xa, ya, 1]
for (xa, ya, _), (xb, yb, _) in zip(x1n, x2n)])
_, _, Vt = np.linalg.svd(A)
F = Vt[-1].reshape(3, 3)
U, S, Vt2 = np.linalg.svd(F) # enforce rank-2 by zeroing the smallest singular value
S[-1] = 0
F = U @ np.diag(S) @ Vt2
F = T2.T @ F @ T1 # undo the normalization
return F / F[2, 2]
F_mine = eight_point_algorithm(x1, x2)
F_cv, _ = cv2.findFundamentalMat(x1, x2, cv2.FM_8POINT)
print('our F:\n', np.round(F_mine, 5))
print('cv2.findFundamentalMat F:\n', np.round(F_cv, 5))
Checking the epipolar constraint directly¶
For true correspondences, $x_2^\top F x_1$ should be (numerically) zero.
def epipolar_residual(F, x1, x2):
x1h = np.hstack([x1, np.ones((len(x1), 1))])
x2h = np.hstack([x2, np.ones((len(x2), 1))])
return np.abs(np.sum(x2h * (F @ x1h.T).T, axis=1))
print(f'mean |x2^T F x1|, our F: {epipolar_residual(F_mine, x1, x2).mean():.2e}')
print(f'mean |x2^T F x1|, cv2 F: {epipolar_residual(F_cv, x1, x2).mean():.2e}')
Visualizing epipolar lines¶
For a handful of points in image 1, we draw their epipolar lines $l_2 = Fx_1$ in image 2, and confirm each true match sits exactly on its line.
w, h = 640, 480
sample_idx = rng.choice(len(x1), 6, replace=False)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(x1[:, 0], x1[:, 1], c='gray', s=15)
axes[0].scatter(x1[sample_idx, 0], x1[sample_idx, 1], c='red', s=40)
axes[0].set_xlim(0, w); axes[0].set_ylim(h, 0)
axes[0].set_title('Image 1: 6 selected points')
axes[1].scatter(x2[:, 0], x2[:, 1], c='gray', s=15)
for i in sample_idx:
a, b, c = F_mine @ np.array([x1[i, 0], x1[i, 1], 1]) # line: a*x + b*y + c = 0
xs = np.array([0, w])
ys = -(a * xs + c) / b
axes[1].plot(xs, ys, linewidth=1)
axes[1].scatter(x2[sample_idx, 0], x2[sample_idx, 1], c='red', s=40, zorder=5, label='true match')
axes[1].set_xlim(0, w); axes[1].set_ylim(h, 0)
axes[1].set_title('Image 2: epipolar lines + true matches')
axes[1].legend(fontsize=8)
plt.tight_layout()
plt.show()
Every true match lands exactly on its predicted line — a point's search in the second image really does collapse from 2D to 1D, even without rectifying the images first.
From fundamental to essential: adding calibration¶
The fundamental matrix works in raw pixel coordinates and folds in the (unknown) camera intrinsics. If we do know the intrinsics $K$ (from camera calibration), we can remove them and work in normalized coordinates, giving the essential matrix:
$$E = K_2^\top F K_1 \qquad \text{(same } K \text{ for both, if it's the same camera)}$$
Unlike $F$ (7 degrees of freedom), $E$ has only 5 — it's built entirely from a relative rotation $R$ and translation direction $t$ between the two cameras: $E = [t]_\times R$. This means $E$ can be decomposed back into $R$ and $t$, which is how you recover camera motion from image correspondences alone.
We estimate $E$ robustly with cv2.findEssentialMat(..., method=cv2.RANSAC, ...), exactly the kind of outlier-rejecting fit built from scratch in Lesson 22.
E_from_F = K.T @ F_mine @ K
E_direct, _ = cv2.findEssentialMat(x1, x2, K, method=cv2.RANSAC, threshold=1.0)
# E is also only defined up to scale -- compare directions, not raw magnitudes
print('E derived from our F (normalized):\n', np.round(E_from_F / np.linalg.norm(E_from_F), 4))
print('E from cv2.findEssentialMat (normalized):\n', np.round(E_direct / np.linalg.norm(E_direct), 4))
Recovering camera motion¶
cv2.recoverPose decomposes $E$ into a rotation and a translation direction (translation magnitude is fundamentally unrecoverable from two uncalibrated-scale views alone — a scene and cameras twice as far apart, moving twice as much, produce identical images. This is the same scale ambiguity familiar from monocular vision).
_, R_estimated, t_estimated, _ = cv2.recoverPose(E_direct, x1, x2, K)
print('true rotation:\n', np.round(R_true, 4))
print('recovered rotation:\n', np.round(R_estimated, 4))
print()
print('true translation direction: ', np.round(t_true / np.linalg.norm(t_true), 4))
print('recovered translation direction: ', np.round(t_estimated.ravel(), 4))
The rotation is recovered essentially exactly, and the translation direction matches up to the expected sign/scale — from feature correspondences alone (Lesson 19), with no other information, we've recovered how the second camera is rotated and (up to scale) positioned relative to the first. This is the starting point for structure-from-motion and visual SLAM.
Exercise¶
- Add pixel noise (e.g. std 0.5) to
x1andx2before running the 8-point algorithm. How much does the mean epipolar residual grow, and does normalizing coordinates (aseight_point_algorithmdoes) actually matter here — try skipping the normalization step and compare. - The epipole in image 2 is the projection of camera 1's center, and satisfies $F e_1 = 0$ for the epipole $e_1$ in image 1 (and $F^\top e_2 = 0$ for the epipole $e_2$ in image 2). Compute both epipoles as the null space of $F$ (via SVD) and check whether they fall inside or outside the visible image region for this camera configuration.
- Increase the rotation angle between the two cameras to 60 degrees and rerun the pose recovery. Does
cv2.recoverPosestill find the correct rotation? At what point would you expect correspondence matching itself (Lesson 19) to become the bottleneck rather than the geometry?