{ "cells": [ { "cell_type": "markdown", "id": "d21c0372", "metadata": {}, "source": [ "# Lesson 28: Epipolar Geometry\n", "\n", "Lesson 22's stereo matching assumed a *rectified* pair, where corresponding points always fall on the same row. In the more general case, stereo cameras may not be rectified. This lesson covers the general two-view relationship between **any** pair of images of a static scene. As we will see, a point's match in the other image isn't just *somewhere* — it's constrained to lie on a particular line. We cover both the calibrated and uncalibrated cases." ] }, { "cell_type": "code", "execution_count": null, "id": "92f2eb83", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "29d2fa5d", "metadata": {}, "source": [ "## The epipolar constraint\n", "\n", "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$ rank-deficient matrix $F$ — the **fundamental matrix** — such that\n", "\n", "$$x_2^\\top F x_1 = 0$$\n", "\n", "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 22's row-restricted search valid — rectification is just the special camera arrangement where every epipolar line happens to be horizontal, which collapses the 1D search to a single row." ] }, { "cell_type": "markdown", "id": "99988f35", "metadata": {}, "source": [ "### Why a line? A geometric picture\n", "\n", "Camera 1 knows only the *direction* to a point, not its depth: the true 3D point $X$ could be anywhere along the ray from $C_1$ through $x_1$. But every point on that ray, together with $C_1$ and $C_2$, lies in a single plane — the **epipolar plane**. Camera 2 sees this plane edge-on, as a single line: wherever $X$ actually sits along the ray, its projection into image 2 always falls on that one line, the epipolar line." ] }, { "cell_type": "code", "execution_count": null, "id": "ede1670d", "metadata": {}, "outputs": [], "source": [ "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", "\n", "C1 = np.array([0.0, 0.0, 0.0])\n", "C2 = np.array([3.0, 0.5, -0.5])\n", "ray_dir = np.array([0.6, 0.8, 3.0])\n", "ray_dir = ray_dir / np.linalg.norm(ray_dir)\n", "depths = [2.0, 3.5, 5.0] # a few candidate depths along camera 1's ray\n", "candidates = [C1 + d * ray_dir for d in depths]\n", "\n", "fig = plt.figure(figsize=(6, 5))\n", "ax = fig.add_subplot(111, projection='3d')\n", "\n", "ax.plot(*zip(C1, C2), c='gray', linestyle='--', linewidth=1)\n", "for X in candidates:\n", " ax.plot(*zip(C1, X), c='tab:blue', linewidth=1)\n", " ax.plot(*zip(C2, X), c='tab:green', linewidth=1)\n", "cand_arr = np.array(candidates)\n", "ax.scatter(cand_arr[:, 0], cand_arr[:, 1], cand_arr[:, 2], c='black', s=30)\n", "ax.text(*candidates[-1], ' possible $X$\\n (unknown depth)', fontsize=8)\n", "\n", "ax.scatter(*C1, c='tab:blue', s=60)\n", "ax.text(*C1, ' $C_1$', fontsize=10)\n", "ax.scatter(*C2, c='tab:green', s=60)\n", "ax.text(*C2, ' $C_2$', fontsize=10)\n", "\n", "plane = Poly3DCollection([[C1, C2, candidates[-1]]], alpha=0.15, facecolor='tab:orange')\n", "ax.add_collection3d(plane)\n", "centroid = (C1 + C2 + candidates[1]) / 3\n", "ax.text(*centroid, 'epipolar\\nplane', fontsize=8, color='tab:orange', ha='center')\n", "\n", "ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')\n", "ax.set_title('Why the match is constrained to a line')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "34a3a47f", "metadata": {}, "source": [ "The three blue segments are all rays from $C_1$ through the same pixel $x_1$, just extended to three different, equally plausible depths — camera 1 alone can't tell them apart. But each candidate, together with $C_1$ and $C_2$, still lies in the one epipolar plane shown, so their projections into camera 2 (the ends of the green segments) all fall along that plane's intersection with image 2: the epipolar line. $F$ is just the algebraic machinery that predicts this line directly from $x_1$, without ever knowing the depth." ] }, { "cell_type": "markdown", "id": "c24278a4", "metadata": {}, "source": [ "## A synthetic two-camera scene\n", "\n", "We build two cameras with known intrinsics $K$ and a known relative rotation/translation, create some random 3D points, project the points into both cameras, and use the resulting correspondences to recover $F$." ] }, { "cell_type": "code", "execution_count": null, "id": "160b7e1f", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(0)\n", "K = np.array([[500, 0, 320], [0, 500, 240], [0, 0, 1]], dtype=np.float64)\n", "\n", "R1, t1 = np.eye(3), np.zeros(3) # camera 1: at the origin, looking down +z\n", "angle = np.radians(15)\n", "R_true = np.array([[np.cos(angle), 0, np.sin(angle)],\n", " [0, 1, 0],\n", " [-np.sin(angle), 0, np.cos(angle)]])\n", "t_true = np.array([0.5, 0.0, 0.1]) # camera 2: rotated 15 deg, shifted along x\n", "\n", "P1 = K @ np.hstack([R1, t1.reshape(3, 1)])\n", "P2 = K @ np.hstack([R_true, t_true.reshape(3, 1)])\n", "\n", "def project(P, points_3d):\n", " homogeneous = np.hstack([points_3d, np.ones((len(points_3d), 1))])\n", " projected = (P @ homogeneous.T).T\n", " return projected[:, :2] / projected[:, 2:3]\n", "\n", "points_3d = rng.uniform(-1, 1, (40, 3)) + np.array([0, 0, 5]) # in front of both cameras\n", "x1 = project(P1, points_3d)\n", "x2 = project(P2, points_3d)\n", "\n", "print(f'{len(points_3d)} 3D points, projected into both cameras')\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))\n", "axes[0].scatter(x1[:, 0], x1[:, 1], c='tab:blue', s=15)\n", "axes[0].set_xlim(0, 640); axes[0].set_ylim(480, 0)\n", "axes[0].set_title('Image 1')\n", "axes[1].scatter(x2[:, 0], x2[:, 1], c='tab:orange', s=15)\n", "axes[1].set_xlim(0, 640); axes[1].set_ylim(480, 0)\n", "axes[1].set_title('Image 2')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "af967249", "metadata": {}, "source": [ "### SVD: Singular value decomposition\n", "\n", "In Lesson 23, we saw that a real symmetric positive-semidefinite matrix $A$ (like a covariance matrix) can be written as $A = P \\Lambda P^\\top$, with $P$ containing the *eigenvectors*, and with $\\Lambda$ diagonal and non-negative containing the *eigenvalues*. The **singular value decomposition (SVD)** generalizes this idea: any matrix $A$ — symmetric or not, square or not — can be written as $A = U \\Sigma V^\\top$, with $U$ and $V$ orthogonal containing the **singular vectors** (left and right, respectively) and $\\Sigma$ diagonal and non-negative containing the **singular values**. By definition, the singular values in $\\Sigma$ are sorted largest to smallest. For our purposes, we focus on two specific facts about the SVD:\n", "\n", "**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 associated with the smallest singular value. This is Lesson 23'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 here is exactly the same \"direction of least fit error\" as the smallest-eigenvalue eigenvector Lesson 23 used.\n", "\n", "**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." ] }, { "cell_type": "markdown", "id": "a6fdf675", "metadata": {}, "source": [ "## The 8-point algorithm\n", "\n", "Each correspondence gives one linear equation in the 9 unknown entries of $F$ (expanding $x_2^\\top F x_1 = 0$). $F$ is only defined up to scale (like the homography in Lesson 24), so it really has just 8 independent unknowns, solvable as a homogeneous least-squares problem via SVD (above) — the same DLT machinery used for the homography fit in Lesson 24. Since $F$ must be exactly rank 2 (a fundamental matrix is always singular), we project the 9-parameter solution back onto the nearest rank-2 matrix, using the SVD trick above.\n", "\n", "(Note that $F$ actually has only 7 true degrees of freedom, because the rank-2 constraint removes one more; however, any 7-point algorithm will be more complicated because it has to take into account a nonlinear constraint; the 8-point algorithm is nice because it is a *linear* algorithm, even though it requires a separate step to restore the rank-2 property.)\n", "\n", "One more wrinkle before we can just hand this to SVD: raw pixel coordinates (in the hundreds) make the linear system badly conditioned numerically. So we first rescale each image's points to be centered at the origin with average distance $\\sqrt2$ from it, solve in these normalized coordinates, then undo the rescaling on the resulting $F$. This conditioning step, shown to make a dramatic difference in practice, is what turns the plain 8-point algorithm into the **normalized 8-point algorithm** (Hartley, 1997) — the version used below, and the one worth reaching for in practice." ] }, { "cell_type": "code", "execution_count": null, "id": "9e269962", "metadata": {}, "outputs": [], "source": [ "def normalize_points(x):\n", " \"\"\"Shift/scale so points are centered at the origin with average distance sqrt(2) -- standard\n", " numerical-conditioning trick (Hartley normalization) for the 8-point algorithm.\"\"\"\n", " mean = x.mean(axis=0)\n", " std = x.std()\n", " T = np.array([[1 / std, 0, -mean[0] / std],\n", " [0, 1 / std, -mean[1] / std],\n", " [0, 0, 1]])\n", " x_h = np.hstack([x, np.ones((len(x), 1))])\n", " return (T @ x_h.T).T, T\n", "\n", "def eight_point_algorithm(x1, x2):\n", " x1n, T1 = normalize_points(x1)\n", " x2n, T2 = normalize_points(x2)\n", "\n", " A = np.array([[xb * xa, xb * ya, xb, yb * xa, yb * ya, yb, xa, ya, 1]\n", " for (xa, ya, _), (xb, yb, _) in zip(x1n, x2n)])\n", " _, _, Vt = np.linalg.svd(A)\n", " F = Vt[-1].reshape(3, 3)\n", "\n", " U, S, Vt2 = np.linalg.svd(F) # enforce rank-2 by zeroing the smallest singular value\n", " S[-1] = 0\n", " F = U @ np.diag(S) @ Vt2\n", "\n", " F = T2.T @ F @ T1 # undo the normalization\n", " return F / F[2, 2]\n", "\n", "F_mine = eight_point_algorithm(x1, x2)\n", "F_cv, _ = cv2.findFundamentalMat(x1, x2, cv2.FM_8POINT)\n", "\n", "print('our F:\\n', np.round(F_mine, 5))\n", "print('cv2.findFundamentalMat F:\\n', np.round(F_cv, 5))" ] }, { "cell_type": "markdown", "id": "40dd50c2", "metadata": {}, "source": [ "### Checking the epipolar constraint directly\n", "\n", "For true correspondences, $x_2^\\top F x_1$ should be (numerically) zero." ] }, { "cell_type": "code", "execution_count": null, "id": "39bbcb90", "metadata": {}, "outputs": [], "source": [ "def epipolar_residual(F, x1, x2):\n", " x1h = np.hstack([x1, np.ones((len(x1), 1))])\n", " x2h = np.hstack([x2, np.ones((len(x2), 1))])\n", " return np.abs(np.sum(x2h * (F @ x1h.T).T, axis=1))\n", "\n", "print(f'mean |x2^T F x1|, our F: {epipolar_residual(F_mine, x1, x2).mean():.2e}')\n", "print(f'mean |x2^T F x1|, cv2 F: {epipolar_residual(F_cv, x1, x2).mean():.2e}')" ] }, { "cell_type": "markdown", "id": "ade675b2", "metadata": {}, "source": [ "### Visualizing epipolar lines\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "ab1802a6", "metadata": {}, "outputs": [], "source": [ "w, h = 640, 480\n", "sample_idx = rng.choice(len(x1), 6, replace=False)\n", "colors = plt.cm.tab10.colors[:len(sample_idx)]\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", "axes[0].scatter(x1[:, 0], x1[:, 1], c='gray', s=15)\n", "axes[0].scatter(x1[sample_idx, 0], x1[sample_idx, 1], c=colors, s=40)\n", "axes[0].set_xlim(0, w); axes[0].set_ylim(h, 0)\n", "axes[0].set_title('Image 1: 6 selected points')\n", "\n", "axes[1].scatter(x2[:, 0], x2[:, 1], c='gray', s=15)\n", "for i, color in zip(sample_idx, colors):\n", " a, b, c = F_mine @ np.array([x1[i, 0], x1[i, 1], 1]) # line: a*x + b*y + c = 0\n", " xs = np.array([0, w])\n", " ys = -(a * xs + c) / b\n", " axes[1].plot(xs, ys, linewidth=1, color=color)\n", "axes[1].scatter(x2[sample_idx, 0], x2[sample_idx, 1], c=colors, s=40, zorder=5,\n", " edgecolors='black', linewidths=0.6, label='true match')\n", "axes[1].set_xlim(0, w); axes[1].set_ylim(h, 0)\n", "axes[1].set_title('Image 2: epipolar lines + true matches')\n", "axes[1].legend(fontsize=8)\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "edb41227", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "3588193c", "metadata": {}, "source": [ "## From fundamental to essential: adding calibration\n", "\n", "The fundamental matrix works in raw pixel coordinates, assuming the camera intrinsics are unknown. If we *do* know the intrinsics $K$ (from camera calibration — Lesson 27), we can remove them and work in normalized coordinates, giving the **essential matrix**:\n", "\n", "$$E = K_2^\\top F K_1 \\qquad \\text{(if both are the same camera, then } K_1 = K_2 = K \\text{)}$$\n", "\n", "Unlike $F$ which has 7 degrees of freedom, $E$ has only 5 degrees of freedom — it's built entirely from a relative rotation $R$ and translation direction $t$ between the two cameras: $E = [t]_\\times R$, where $[t]_\\times$ is the $3\\times3$ skew-symmetric \"cross-product matrix\" of $t$, defined so that $[t]_\\times v = t \\times v$ for any vector $v$. This means $E$ can be *decomposed* back into $R$ and $t$, which is how you recover camera motion from image correspondences alone.\n", "\n", "We estimate $E$ robustly with `cv2.findEssentialMat(..., method=cv2.RANSAC, ...)`, exactly the kind of outlier-rejecting fit built from scratch in Lesson 23." ] }, { "cell_type": "code", "execution_count": null, "id": "05f38509", "metadata": {}, "outputs": [], "source": [ "E_from_F = K.T @ F_mine @ K\n", "E_direct, _ = cv2.findEssentialMat(x1, x2, K, method=cv2.RANSAC, threshold=1.0)\n", "\n", "# E is also only defined up to scale -- compare directions, not raw magnitudes\n", "print('E derived from our F (normalized):\\n', np.round(E_from_F / np.linalg.norm(E_from_F), 4))\n", "print('E from cv2.findEssentialMat (normalized):\\n', np.round(E_direct / np.linalg.norm(E_direct), 4))" ] }, { "cell_type": "markdown", "id": "7a595543", "metadata": {}, "source": [ "### Recovering camera motion\n", "\n", "`cv2.recoverPose` decomposes $E$ into a rotation and a translation *direction* — not a full translation vector, since the magnitude is fundamentally unrecoverable from two views alone. A scene twice as large, viewed by cameras with twice the baseline between them, produces identical images; this is the same scale ambiguity familiar from monocular vision." ] }, { "cell_type": "code", "execution_count": null, "id": "5b25c82c", "metadata": {}, "outputs": [], "source": [ "_, R_estimated, t_estimated, _ = cv2.recoverPose(E_direct, x1, x2, K)\n", "\n", "print('true rotation:\\n', np.round(R_true, 4))\n", "print('recovered rotation:\\n', np.round(R_estimated, 4))\n", "print()\n", "print('true translation direction: ', np.round(t_true / np.linalg.norm(t_true), 4))\n", "print('recovered translation direction: ', np.round(t_estimated.ravel(), 4))" ] }, { "cell_type": "markdown", "id": "2dae4088", "metadata": {}, "source": [ "The rotation is recovered, and the translation direction matches up to the expected sign/scale. From feature correspondences alone (Lesson 20), we've recovered the relative rotation and translation (up to scale) between the two cameras. This is the starting point for structure-from-motion (Lesson 29)." ] }, { "cell_type": "markdown", "id": "b5962a05", "metadata": {}, "source": [ "### Exercises\n", "\n", "1. Add pixel noise (e.g. std 0.5) to `x1` and `x2` before running the 8-point algorithm. How much does the mean epipolar residual grow, and does normalizing coordinates (as `eight_point_algorithm` does) actually matter here — try skipping the normalization step and compare.\n", "2. 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.\n", "3. Increase the rotation angle between the two cameras to 60 degrees and rerun the pose recovery. Does `cv2.recoverPose` still find the correct rotation? At what point would you expect correspondence matching itself (Lesson 20) to become the bottleneck rather than the geometry?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }