From Pixels to 3D: Features, RANSAC & Structure from Motion

Chapter 1: An Image Is a Function

Before we can find corners or match points, we need the right mental model of an image. To a robot, a grayscale image is not a picture — it is a function. Stanford's Lecture 9 writes it precisely:

Read that aloud: I takes a position (u, v) — a column and a row in the image — and returns a single number, the intensity (brightness) at that spot, between 0 (black) and L (white, usually 255). A color image is just three such functions stacked, one each for red, green, blue. From here on we work in grayscale; everything generalizes.

Because the function lives on a discrete grid of pixels, we store it as a matrix: row index v, column index u, and the entry I(u, v) is the brightness.

The gradient: how brightness changes

The image has no formula, so we cannot differentiate it the calculus-class way. Instead we approximate the derivative with a finite difference — literally subtracting neighboring pixels. The centered difference from Lecture 9 is:

Here I_u is the horizontal gradient (how fast brightness rises as you walk right) and I_v is the vertical gradient. Stack them and you get the gradient vector at each pixel:

This little vector points uphill in brightness, in the direction of fastest increase, and its length tells you how steep the change is. A bright edge has a long gradient vector pointing across the edge; a flat region has a near-zero gradient.

Worked example: gradient of a 3×3 patch by hand

Let's compute one gradient with real numbers. Suppose a tiny patch of brightness looks like this (a vertical dark-to-light edge), and we want the gradient at the center pixel (the 90):

Apply the centered differences at the center pixel:

So ∇I = (50, 30). The gradient magnitude — the steepness, used for edge detection — is its length:

And the gradient direction (which way is brightness climbing fastest) is:

The brightness climbs mostly to the right and a bit downward. An edge runs perpendicular to the gradient, so this edge is roughly vertical — exactly what the dark-left / light-right pattern shows. We just detected an edge with three subtractions.

From scratch in Python

python
import numpy as np

def gradients(I):
    # I: float grayscale image, shape (H, W)
    Iu = np.zeros_like(I)
    Iv = np.zeros_like(I)
    # centered difference; skip the 1-pixel border
    Iu[:, 1:-1] = (I[:, 2:] - I[:, :-2]) / 2.0   # horizontal
    Iv[1:-1, :] = (I[2:, :] - I[:-2, :]) / 2.0   # vertical
    mag = np.sqrt(Iu**2 + Iv**2)              # edge strength
    ang = np.arctan2(Iv, Iu)                  # edge orientation
    return Iu, Iv, mag, ang

And the idiomatic, fast version uses a Sobel filter (a smoothed derivative) from OpenCV:

python
import cv2
Iu = cv2.Sobel(I, cv2.CV_64F, 1, 0, ksize=3)   # d/du, blurred
Iv = cv2.Sobel(I, cv2.CV_64F, 0, 1, ksize=3)   # d/dv, blurred
mag = cv2.magnitude(Iu, Iv)

Run the gradient field below. Move the slider to change the edge angle in a synthetic image; the little arrows show ∇I at each cell — long across the edge, short along it.

Chapter 2: Finding Corners — The Harris Detector

Edges are useful, but for matching across views they have a fatal flaw. Slide a window along an edge and the picture barely changes — one stretch of edge looks like any other stretch. You cannot tell where on the edge you are. That ambiguity is the aperture problem, and it makes edges poor landmarks.

What you want is a corner: a point where two edges meet, so that brightness changes in every direction. Slide a window even slightly off a corner and the picture changes a lot, in a unique way. That is what makes corners repeatable (findable in both views) and distinctive (matchable) — the two properties Lecture 9 demands of a good feature.

The intuition: shift a window and watch the change

Stanford's geometric picture: take a small window, shift it by a tiny amount (Δu, Δv), and measure how much the patch's intensity changes (sum of squared differences). Three cases:

• Flat region: shifting any direction changes nothing. No feature.
• Edge: shifting along the edge changes nothing; shifting across it changes a lot. One direction only.
• Corner: shifting any direction changes the patch a lot. This is the winner.

The change for a small shift turns out (after a Taylor expansion) to depend on a single 2×2 matrix built from the gradients in the window — the structure tensor (also called the second-moment matrix), M:

Each entry is a sum over all pixels in the window. Define the three sums:

This matrix is a compact summary of which directions have brightness change inside the window. Its two eigenvalues λ₁, λ₂ measure the change along the two principal directions:

• Both eigenvalues small → flat (no change anywhere).
• One large, one small → edge (change in one direction).
• Both large → corner (change in every direction).

The cornerness score — no eigenvalues needed

Computing eigenvalues per pixel is expensive. Harris's trick: a score that is large only when both eigenvalues are large, built from cheap quantities — the determinant and trace of M:

where det(M) = λ₁λ₂ = AB − C², tr(M) = λ₁ + λ₂ = A + B, and k is a small constant, typically 0.04–0.06. Why does this work? The determinant is the product of eigenvalues — it is only large when both are large. The trace penalty subtracts off cases where one eigenvalue dominates (an edge), where the trace squared blows up but the product does not. So:

• R > 0 and large → corner.
• R < 0 → edge (one eigenvalue ≫ the other).
• |R| small → flat.

Worked example: cornerness of a tiny patch by hand

Take a window with three pixels whose gradients we computed (mix of horizontal and vertical change — a corner-like patch). Gradients (I_u, I_v) at the three pixels are (4, 0), (0, 3), (2, 2). Build the three sums:

R = 200.4 > 0 and large → this is a corner. Now contrast with a pure-edge patch where all gradients are horizontal: (4, 0), (5, 0), (3, 0). Then A = 16+25+9 = 50, B = 0, C = 0:

Negative → correctly rejected as an edge. The determinant collapsed to zero because all change was in one direction. Harris just told an edge from a corner with grade-school arithmetic.

From scratch in Python

python
import numpy as np
from scipy.ndimage import uniform_filter

def harris(I, k=0.04, win=3):
    Iu = np.zeros_like(I); Iv = np.zeros_like(I)
    Iu[:, 1:-1] = (I[:, 2:] - I[:, :-2]) / 2.0
    Iv[1:-1, :] = (I[2:, :] - I[:-2, :]) / 2.0
    # sum the products over a window == the structure tensor entries A,B,C
    A = uniform_filter(Iu*Iu, win)
    B = uniform_filter(Iv*Iv, win)
    C = uniform_filter(Iu*Iv, win)
    detM = A*B - C*C
    trM  = A + B
    R = detM - k*trM**2                # cornerness per pixel
    return R

# corners = local maxima of R above a threshold (non-max suppression)

And the one-liner most engineers actually call:

python
import cv2
R = cv2.cornerHarris(I.astype('float32'), blockSize=2, ksize=3, k=0.04)
corners = np.argwhere(R > 0.01*R.max())   # threshold + (you still need NMS)

Below, the cornerness heatmap is computed live on a synthetic scene with a square. Bright spots are high R. Watch the four corners of the square light up red while the edges stay dark — exactly what the determinant predicts. Drag k to see how it trades corners against edges.

Chapter 3: Describing & Matching Features

Harris hands us a list of corner locations: "interesting point at pixel (412, 88)." But a location alone cannot be matched across two photos. To say "this corner in image A is that corner in image B," each corner needs a fingerprint — a vector of numbers that summarizes the patch around it. That fingerprint is a descriptor.

Lecture 9 lays out the design requirements. A good descriptor should be:

• Invariant — the same physical point should produce (nearly) the same descriptor despite changes in viewpoint, scale, rotation, and lighting.
• Distinctive — different points should produce different descriptors, or you cannot tell them apart.

The naive descriptor and why it fails

The simplest descriptor: grab the n×m window of raw pixel intensities around the corner and flatten it into a vector. Normalize it (subtract the mean, divide by the standard deviation) and it becomes invariant to linear brightness changes. Stanford lists its three fatal drawbacks:

1. Sensitive to pose. Rotate the camera a few degrees and every pixel shifts; the raw window changes completely.
2. Sensitive to scale. Move closer and the patch covers different physical area.
3. Poorly distinctive. Two different smooth patches look alike pixel-for-pixel.

The SIFT idea: a histogram of gradient orientations

The SIFT descriptor (Scale-Invariant Feature Transform) is the workhorse. Its key idea, in one breath: describe a keypoint by the pattern of gradient orientations in its neighborhood, not the raw intensities. Orientations survive lighting changes (a brighter image has the same gradient directions, just longer arrows), and binning them into a histogram tolerates small shifts.

The construction, simplified:

1. Take a patch (say 16×16) around the keypoint and compute the gradient (magnitude, orientation) at every pixel.
2. Assign a dominant orientation to the whole patch (the peak of its orientation histogram) and rotate the patch to that angle. Now the descriptor is rotation-invariant: a corner described the same way no matter how the camera was tilted.
3. Split the patch into a 4×4 grid of cells. In each cell, build an 8-bin histogram of gradient orientations, weighted by magnitude.
4. Concatenate: 4×4 cells × 8 bins = 128 numbers. Normalize to unit length for lighting invariance. That 128-dimensional vector is the SIFT descriptor.

The genius is in step 3: a histogram throws away exact pixel positions within each cell, keeping only "how much gradient points which way." That tolerance is precisely what makes it robust to the small geometric wobble between two views.

Matching: nearest descriptor + the Lowe ratio test

Now we have a set of 128-D descriptors in image A and another set in image B. To match a descriptor d_A, find the nearest descriptor in B by Euclidean distance:

But "nearest" is not enough. Smooth regions and repetitive textures produce many near-identical descriptors, and the closest one is often a false match. David Lowe's fix — the single most important line in feature matching — is the ratio test: find the two nearest neighbors in B (distances d₁ ≤ d₂) and keep the match only if the best is clearly better than the runner-up:

Worked example: the ratio test by hand

Descriptor d_A's two nearest neighbors in image B are at distances d₁ = 0.18 and d₂ = 0.21. The ratio is:

The best match is barely better than the second-best — ambiguous, probably a repetitive texture, so we throw it away. Now another descriptor with d₁ = 0.12 and d₂ = 0.40:

The best match is far better than any alternative — a confident, unique correspondence. The ratio test is a cheap, brutally effective outlier filter that discards most bad matches before geometry ever sees them. But it does not catch all of them — which is why we still need RANSAC (Chapter 6).

From scratch in Python

python
import numpy as np

def match_ratio(descA, descB, tau=0.8):
    matches = []
    for i, da in enumerate(descA):
        # Euclidean distance from da to every descriptor in B
        dists = np.linalg.norm(descB - da, axis=1)
        order = np.argsort(dists)            # nearest first
        d1, d2 = dists[order[0]], dists[order[1]]
        if d1 < tau * d2:                     # Lowe ratio test
            matches.append((i, order[0]))    # (index in A, index in B)
    return matches

And the library version with OpenCV's SIFT + brute-force matcher:

python
import cv2
sift = cv2.SIFT_create()
kpA, descA = sift.detectAndCompute(imgA, None)
kpB, descB = sift.detectAndCompute(imgB, None)
bf = cv2.BFMatcher()
pairs = bf.knnMatch(descA, descB, k=2)
good = [m for m, n in pairs if m.distance < 0.8*n.distance]   # ratio test

Try the matcher below. Two views of a scene with several corners. Slide the ratio threshold τ: low τ keeps only the most confident matches (few lines, all correct); high τ keeps more, but lets a few wrong (red) lines slip through.

Chapter 4: Epipolar Geometry — The Geometry of Two Views

Matching by descriptors is a 2D search over millions of pixels — slow and error-prone. But two views are not independent: they are rigidly tied by geometry. Once you know how the two cameras are placed relative to each other, the match for a point in image A is forced to lie on a single line in image B. That constraint, the crown jewel of two-view vision, is epipolar geometry.

The setup

Two cameras with centers O and O′ both look at a 3D point P. It projects to pixel p in image 1 and p′ in image 2. Now stare at the five points P, O, O′, p, p′. They are coplanar — they all lie on one plane, the epipolar plane. Why? Because p is on the ray Op, p′ is on the ray O′p′, and both rays pass through P — so the two rays plus the line joining the camera centers (the baseline OO′) all live in the plane they define.

Some vocabulary you will see everywhere:

• The baseline is the line connecting the two camera centers O and O′.
• The epipole in each image is where the baseline pierces that image plane — equivalently, the image of the other camera's center. Call them e and e′.
• An epipolar line is where the epipolar plane meets an image plane. Every epipolar line in an image passes through that image's epipole.

The fundamental matrix F

This coplanarity condition can be written as one beautifully compact equation. With p and p′ written in homogeneous coordinates (a pixel (u, v) becomes the 3-vector (u, v, 1)), the epipolar constraint is:

Here F is the fundamental matrix, a 3×3 matrix that encodes everything about the geometric relationship between the two views — the relative rotation, the relative translation, and both cameras' internal calibration, all baked into nine numbers. Let's unpack the equation:

• p is a point in image 1 (homogeneous 3-vector).
• F p is a 3-vector that represents the epipolar line l′ in image 2 (a line au + bv + c = 0 is stored as the vector (a, b, c)).
• p′^⊤(Fp) = 0 says: the matching point p′ lies on that line. (A point lies on a line exactly when their dot product is zero.)

Fundamental F vs. essential E

There are two cousins, and the difference is whether you have calibrated cameras:

The exact relationship, from Lecture 18 (with K the camera intrinsic matrix, R the relative rotation, and [t]_× the skew-symmetric matrix encoding the translation as a cross product):

The essential matrix E is the prize: once you have it (and you have K), you can decompose it to recover the relative rotation R and translation t between the two cameras — the camera motion. That decomposition is the engine of visual odometry (Chapter 9).

Click a point in the left image below. The constraint draws the matching epipolar line in the right image — the only place its match can live. Drag the second camera (baseline slider) to watch the epipoles and the family of lines pivot.

Chapter 5: The 8-Point Algorithm — Estimating F

The fundamental matrix is so useful that the obvious question is: how do we find it from a pair of images? We don't know the camera motion in advance. But we do have something better — a pile of corresponding points (p_i, p′_i) from Chapter 3. Each correspondence gives us one equation about F, and enough of them pin it down. This is the 8-point algorithm.

Turning the constraint into a linear equation

Write the constraint p′^⊤ F p = 0 for one correspondence, with p = (u, v, 1) and p′ = (u′, v′, 1) and the unknown matrix written out:

Multiplying it all out, the single scalar equation p′^⊤Fp = 0 becomes linear in the nine unknowns:

Stack the nine unknowns into a vector f = (f₁₁, …, f₃₃). Then each correspondence contributes one row a_i of known coefficients:

With n correspondences we stack the rows into an n×9 matrix A and demand:

Why "8 points"?

F has nine entries, but it is only defined up to scale (multiplying every entry by 2 leaves the constraint p′^⊤Fp = 0 unchanged). So there are really only 8 free parameters, and you need 8 correspondences to fix them. With exactly 8 points you solve the system directly; with more, you solve a least-squares problem (overdetermined, more robust).

Solving Af = 0 by least squares (the SVD)

We cannot just solve Af = 0 the usual way, because f = 0 is a trivial useless answer. We want the nonzero f (of unit length) that makes Af as close to zero as possible:

The answer is a classic linear-algebra fact: it is the singular vector of A with the smallest singular value — the last column of V in the SVD A = UΣV^⊤. Reshape that 9-vector back into a 3×3 matrix and you have an estimate F̂.

Two crucial fix-ups

The raw estimate F̂ needs two corrections, both straight from Lecture 18:

1. Normalization (Hartley). Pixel coordinates have wildly different scales (e.g. u′u can be hundreds of thousands while the last entry is 1). This makes A badly conditioned and the answer noise-sensitive. The fix: before building A, translate and scale each image's points so they are centered at the origin with average distance √2. Solve, then undo the transform. This is not optional — it is the difference between garbage and a usable F.
2. Rank-2 enforcement. A true fundamental matrix is singular (det F = 0, rank 2), but noisy data gives a full-rank estimate. Force it: take the SVD of F̂ = UΣV^⊤, set the smallest singular value in Σ to zero, and rebuild. This is the closest rank-2 matrix (in Frobenius norm) to your estimate.

Worked setup: building one row of A

Suppose one correspondence is p = (100, 50, 1) in image 1 and p′ = (110, 48, 1) in image 2. Its row of A is:

Notice the disastrous spread: 11000 down to 1. That is exactly why Hartley normalization — rescaling the points first — is essential. Eight such rows stack into A, and the smallest-singular-vector gives f.

From scratch in Python

python
import numpy as np

def eight_point(pts1, pts2):
    # pts1, pts2: arrays of shape (n, 2), n >= 8 correspondences
    n = pts1.shape[0]
    A = np.zeros((n, 9))
    for i in range(n):
        u, v = pts1[i]; up, vp = pts2[i]
        A[i] = [up*u, up*v, up, vp*u, vp*v, vp, u, v, 1]
    # smallest singular vector of A -> f
    _, _, Vt = np.linalg.svd(A)
    F = Vt[-1].reshape(3, 3)
    # enforce rank 2: zero the smallest singular value
    U, S, Vt2 = np.linalg.svd(F)
    S[2] = 0
    F = U @ np.diag(S) @ Vt2
    return F / F[2, 2]   # normalize scale

(A real implementation wraps this in Hartley normalization.) The library version:

python
import cv2
F, mask = cv2.findFundamentalMat(pts1, pts2, cv2.FM_RANSAC)   # normalization + rank-2 + RANSAC

Notice that flag: FM_RANSAC. Real point lists are full of bad matches that no ratio test caught, and a single outlier destroys a least-squares fit. The cure is the algorithm in the next chapter.

Chapter 6: RANSAC — Fitting Through a Storm of Outliers

Here is the crisis. Your matcher returns 200 correspondences. Maybe 140 are correct (inliers); the other 60 are nonsense (outliers) — repetitive textures, glare, mismatched windows that survived the ratio test. Now fit a fundamental matrix to all 200 with least squares. The result is ruined. Least squares minimizes squared error, so a few wild outliers, squared, dominate the fit. One bad match can drag the answer off a cliff.

We need an estimator that ignores the bad data instead of averaging it in. That estimator is RANSAC — RANdom SAmple Consensus. It is one of the most important algorithms in all of robotics and vision, and the idea is shockingly simple.

The algorithm, step by step

To fit a model needing a minimal sample of s points (e.g. a line needs s = 2, a fundamental matrix needs s = 8):

Then, as a final polish, refit the best model using all its inliers by least squares — now safe, because the outliers have been removed from the pool. A point is an inlier if its error from the model is below a threshold ε (for an epipolar fit, the distance from the matched point to its predicted epipolar line).

How many iterations? The N formula

How many random samples N do we need to be confident? We want at least one sample that is all inliers. Let w be the fraction of inliers in the data. Then:

• The probability one random point is an inlier is w.
• The probability all s points in a sample are inliers is w^s.
• So the probability a sample is "bad" (has at least one outlier) is 1 − w^s.
• The probability that all N samples are bad is (1 − w^s)^N.

We want that failure probability to be small, equal to 1 − p where p is our desired success confidence (e.g. 0.99). Set (1 − w^s)^N = 1 − p and take logs:

Worked example: N for line fitting with 50% outliers

Fit a line (minimal sample s = 2), inlier fraction w = 0.5 (half the data is junk), desired confidence p = 0.99. Plug in:

Stunning. Even with half the data being garbage, just 16 random two-point samples give 99% confidence of hitting one clean sample. That is why RANSAC is everywhere — it is cheap and almost magically robust to outliers.

The fundamental matrix, RANSAC-ified

Combining Chapters 5 and 6 gives the production recipe used in every visual-SLAM front end:

1. Randomly pick 8 correspondences.
2. Fit F with the 8-point algorithm.
3. Count inliers: correspondences where p′ lies within ε pixels of its predicted epipolar line Fp.
4. Repeat N times; keep the F with the most inliers.
5. Refit F using all inliers; discard the outliers forever.

From scratch in Python (line fit)

python
import numpy as np

def ransac_line(pts, eps=1.0, p=0.99, w=0.5):
    s = 2
    N = int(np.ceil(np.log(1-p) / np.log(1 - w**s)))   # iteration budget
    best_inliers, best_model = [], None
    for _ in range(N):
        i, j = np.random.choice(len(pts), 2, replace=False)   # 1. sample
        (x1,y1),(x2,y2) = pts[i], pts[j]
        a, b, c = y2-y1, x1-x2, x2*y1 - x1*y2            # 2. fit line ax+by+c=0
        norm = np.hypot(a, b) + 1e-9
        dist = np.abs(a*pts[:,0] + b*pts[:,1] + c) / norm   # 3. point-to-line dist
        inliers = np.where(dist < eps)[0]
        if len(inliers) > len(best_inliers):              # 4. keep best
            best_inliers, best_model = inliers, (a,b,c)
    # 5. refit on all inliers (least squares) -- omitted for brevity
    return best_model, best_inliers

Watch it run. Points are mostly on a line plus a cloud of outliers. RANSAC repeatedly samples two points (the thin candidate line), counts inliers, and locks onto the consensus line in red — while plain least squares (gray) gets yanked toward the outliers. Crank the outlier slider and watch least squares fail while RANSAC holds.

Chapter 7: Triangulation — From Two Rays to a 3D Point

We have come a long way: detected corners, matched them, estimated the geometry with RANSAC, and (by decomposing the essential matrix) recovered the relative pose (R, t) between the two cameras. Now we cash in. Given a matched pair (p, p′) and the camera poses, triangulation recovers the 3D point P.

The clean idea: intersect two rays

Each pixel back-projects to a ray. Pixel p defines the ray Op from camera 1's center; pixel p′ defines the ray O′p′ from camera 2's center. The 3D point lives where the two rays cross. In a perfect world they intersect exactly — that is the geometry of Chapter 0 made precise.

The catch: noisy rays don't quite meet

In reality, pixels are measured with noise, so the two rays are skew — they pass close but never exactly touch. Lecture 18's fix: instead of demanding an exact intersection, find the 3D point Q whose projections q, q′ into the two images are as close as possible to the measured p, p′. We minimize the re-projection error:

where q(Q) is where the candidate 3D point Q would land in image 1. Re-projection error — "how far off is my reconstruction when I project it back into the photos" — is the universal currency of multi-view geometry. We will see it again in bundle adjustment.

The linear (DLT) method

A point projects via its camera matrix M (a 3×4 matrix combining intrinsics and pose) as p ∼ M Q, where ∼ means "equal up to scale" and Q = (X, Y, Z, 1) is the homogeneous 3D point. The "up to scale" means the cross product p × MQ = 0. Each image gives two independent linear equations in the unknown Q; two images give four equations, and — just like the 8-point algorithm — we solve the resulting homogeneous system with an SVD (smallest singular vector). This is the Direct Linear Transform (DLT).

Worked example: triangulating a point with rectified cameras

To see real numbers without 3×4 matrix grind, take the simplest case — two parallel cameras (the stereo setup of Chapter 8), focal length f = 500 px, baseline B = 0.10 m, with the left camera at the origin. A point projects to u_L = 320 in the left image and u_R = 300 in the right. The disparity is:

Depth (derived in full next chapter) is Z = fB/d:

Now recover X and Y by back-projecting the left pixel (with principal point at, say, c_u = 320, c_v = 240, and v_L = 240):

So P = (0, 0, 2.5) m — straight ahead, 2.5 meters out. We just turned two pixel measurements into a 3D point with arithmetic.

The scale ambiguity

One profound subtlety from Lecture 18: with a single camera moving (structure from motion), you can recover the scene shape and the camera motion only up to an unknown overall scale. The reason: the essential matrix gives the translation direction t but not its length. Double the whole scene and double the camera baseline, and the images are pixel-for-pixel identical. You cannot tell a dollhouse photographed up close from a real house photographed far away.

From scratch (linear triangulation)

python
import numpy as np

def triangulate(p1, p2, M1, M2):
    # p1, p2: pixel (u, v) in each image; M1, M2: 3x4 camera matrices
    u1, v1 = p1; u2, v2 = p2
    A = np.array([
        u1*M1[2] - M1[0],
        v1*M1[2] - M1[1],
        u2*M2[2] - M2[0],
        v2*M2[2] - M2[1],
    ])                                  # 4x4 from p x MQ = 0
    _, _, Vt = np.linalg.svd(A)
    Q = Vt[-1]                          # smallest singular vector
    return Q[:3] / Q[3]                  # de-homogenize -> (X, Y, Z)

And the library one-liner:

python
import cv2
pts4d = cv2.triangulatePoints(M1, M2, pts1.T, pts2.T)   # 4xN homogeneous
pts3d = (pts4d[:3] / pts4d[3]).T                       # Nx3

Drag the matched point in the demo. Two camera rays redraw and intersect; the 3D depth readout updates. Add measurement noise to see the rays go skew and the reconstruction wobble — the re-projection-error minimizer picks the midpoint.

Chapter 8: Stereo Vision — Depth From a Rigid Pair

Triangulation is easiest when the two cameras are arranged in the friendliest possible way: pointing the same direction, offset horizontally, like your two eyes. This is a stereo rig, and it is so common (every depth camera, every self-driving stereo pair) that it earns its own clean formula.

Rectification: make the epipolar lines horizontal

If the two image planes are exactly parallel and the cameras differ only by a horizontal shift, something wonderful happens: every epipolar line becomes a horizontal row, and corresponding points have the same v-coordinate (same row) in both images. The 1D epipolar search of Chapter 4 becomes a search along a single image row — trivial.

Real rigs are never perfectly aligned, so we rectify: apply a projective warp to both images that simulates a perfectly parallel pair. Lecture 18: "from now on, we assume rectified image pairs." After rectification, matching point p′ for p sits on the same row, shifted left by some amount.

Disparity and the depth formula

That horizontal shift is the disparity:

where u_L and u_R are the column of the same point in the left and right images. Now the payoff. By similar triangles (the triangle from the 3D point to the two camera centers is similar to the triangle from the point to the two image projections), the depth is:

Three letters, each earning its place:

• f = focal length in pixels (from calibration, Lesson 8).
• B = baseline, the physical distance between the two cameras (measured once, fixed).
• d = disparity in pixels (measured per point).

Worked example: depth from disparity, with numbers

A stereo camera has focal length f = 700 px and baseline B = 0.12 m. A point appears at column u_L = 410 in the left image and u_R = 380 in the right. Step by step:

The object is 2.8 meters away. Now move it closer so the disparity doubles to d = 60:

Double the disparity, halve the depth — the inverse relationship in action. And a far object with d = 3 px:

Notice the danger: at large depth, a tiny disparity error (3 vs 4 px) swings Z from 28 m to 21 m — a 7-meter jump from one pixel. Stereo is precise up close and increasingly uncertain far away.

The baseline trade-off

This is the same tension as Chapter 0's baseline slider: wider baseline = more depth information but harder correspondence. Stereo rig designers tune B to the working range.

From scratch in Python

python
import numpy as np

def depth_from_disparity(disparity, f, B):
    # disparity: HxW pixel-shift map; f: focal px; B: baseline m
    Z = np.where(disparity > 0, f * B / disparity, np.inf)
    return Z          # per-pixel depth in meters

And computing the disparity map with OpenCV's block matcher:

python
import cv2
stereo = cv2.StereoBM_create(numDisparities=64, blockSize=15)
disp = stereo.compute(left_gray, right_gray).astype(np.float32) / 16.0
Z = f * B / disp           # meters

Drag the object in the rectified stereo demo. Both views update, the disparity bar and depth Z = fB/d recompute live. Watch the disparity shrink as you push the object away.

Chapter 9: Structure from Motion — The Full Pipeline

Stereo had it easy: a fixed rig with a known baseline. Structure from Motion (SfM) is the harder, more general problem from Lecture 18: given m images of n fixed 3D points, taken from unknown positions (maybe even by different cameras at different times), recover both:

• the m camera poses (the motion), and
• the n 3D points (the structure).

This is the algorithm behind a single moving camera reconstructing its surroundings — the core of monocular visual odometry and the front-end of visual SLAM. Every chapter you have read is a stage in it.

The two-view SfM pipeline, assembled

Lecture 18's algebraic two-view recipe, now that you own every piece:

Walk through it once in words. Detect corners in image A and B. Describe and match them, ratio-testing out the obvious junk. Feed the surviving matches to RANSAC + 8-point to get a clean fundamental matrix and an inlier set. Convert F to the essential matrix using the calibration K, and decompose it into the relative rotation and translation. Now both cameras live in one frame, so triangulate every inlier match into a 3D point. You have structure and motion — up to the scale ambiguity of Chapter 7.

Bundle adjustment: the final polish

The two-view estimate is good but not optimal — errors accumulate, and with more than two views the camera poses and 3D points must all agree simultaneously. Bundle adjustment is the global refinement: jointly nudge all camera poses and all 3D points to minimize the total re-projection error across every image:

Read it: for every point j seen in every image i, project the current 3D estimate into that image and compare to the actually-observed pixel q_ij. Sum all those squared errors and minimize over everything at once. It is a large nonlinear least-squares problem, solved with Levenberg-Marquardt. The name evokes "bundles" of light rays converging on each camera, adjusted until they best agree.

Scale drift — the monocular curse

String many two-view steps together (frame 1→2, 2→3, 3→4, …) to track a camera over a long trajectory, and each step's unknown scale (Ch 7) is estimated slightly differently. The errors compound: the reconstructed trajectory slowly grows or shrinks relative to reality. This is scale drift, and it is the defining weakness of monocular visual odometry.

SfM as visual odometry

Lecture 18 closes on the robotics payoff: visual odometry — estimating the robot's motion from its camera stream. A monocular camera gives motion up to scale; a stereo camera gives motion in absolute scale directly (every frame is metric). This is how a robot with cameras knows it moved 1.2 m forward and turned 15° — by running this exact detect→match→estimate→triangulate loop on consecutive frames. Features and SfM feed straight into Perception, and their pose estimates feed Localization & SLAM (Lesson 11).

python
# Two-view SfM, end to end (OpenCV)
import cv2, numpy as np
# 1-2. detect + describe + match (ratio-tested)  --> pts1, pts2  (Ch 2-3)
# 3.   robust essential matrix (RANSAC built in)
E, mask = cv2.findEssentialMat(pts1, pts2, K, method=cv2.RANSAC, prob=0.999)
# 4.   recover relative pose (R, t) by decomposing E
_, R, t, mask = cv2.recoverPose(E, pts1, pts2, K)
# 5.   build camera matrices and triangulate the inliers
M1 = K @ np.hstack([np.eye(3), np.zeros((3,1))])
M2 = K @ np.hstack([R, t])
pts4d = cv2.triangulatePoints(M1, M2, pts1.T, pts2.T)
cloud = (pts4d[:3] / pts4d[3]).T          # 3D point cloud (up to scale)
# 6.   bundle adjustment over many frames -> refined poses + points

Chapter 10: Showcase — RANSAC Fights the Outliers, Live

This is the payoff. One interactive sandbox that ties together the two ideas at the heart of robust geometry: least squares collapses under outliers, and RANSAC does not. You control the outlier fraction; you watch both estimators react in real time; and the iteration budget N = log(1−p)/log(1−w^s) recomputes as you drag.

The scene: a cloud of 2D points that mostly lie on a hidden true line (the inliers), plus a tunable fraction of pure-noise points scattered anywhere (the outliers). Two estimators race to recover the line:

• Least squares (gray): fits all points by minimizing squared vertical error. Every outlier, squared, yanks it.
• RANSAC (red): samples random pairs, scores by consensus, keeps the best, refits on inliers.

Press Run to animate RANSAC sampling: each frame draws the current random 2-point candidate (cyan) and highlights its inliers. The best-so-far consensus line firms up in red. The status panel reports the live inlier count, the consensus count, the running best, and the theoretical iteration budget.

This sandbox is, in miniature, what a robot does when it estimates its motion from a camera: its matched features are full of outliers, and RANSAC is what lets it find the true geometry anyway. Swap "line" for "fundamental matrix" and "point-to-line distance" for "distance to epipolar line," and you have the actual SfM front end.

Chapter 11: Connections & Cheat-Sheet

You started with a flat grid of pixels and ended with a 3D reconstruction and a recovered camera motion. Let's consolidate, place this lesson in the autonomy stack, and point to where it leads.

Where this sits on the stack

Features and SfM are squarely in Perception — they turn raw pixels (Lesson 8's camera images) into geometric meaning: 3D points and relative poses. Their outputs flow upward into Localization & SLAM (Lesson 11), where the camera poses become the robot's trajectory and the 3D points become a feature map. The feature/landmark map from Lesson 1's map zoo? This is the lesson that produces it.

Cheat-sheet

Related lessons on Engineermaxxing

These go deeper on what we touched, and on what comes next:

• AA274A Lesson 8 — Camera Models & Calibration — the pinhole model, intrinsics K, and why one image gives only rays. The direct prerequisite to this lesson.
• AA274A Lesson 11 — SLAM & Factor Graphs — where these feature matches and poses become a full map with loop closure that finally tames scale drift.
• Multi-View Geometry — a deeper dive on epipolar geometry, the fundamental/essential matrices, and homographies.
• Optical Flow — tracking how every pixel moves between frames, the dense cousin of sparse feature matching.
• Visual Navigation — Two-View Geometry — the same epipolar machinery framed for aerial/visual navigation.
• Visual Navigation — RANSAC — RANSAC in full, with adaptive iteration counts and model variants.