Lesson 23: Projective Geometry¶

Lesson 8 built a hierarchy of 2D transforms — Euclidean, similarity, affine — and noted that all three preserve parallel lines. This lesson covers the next, most general step: the projective transform (homography), which models what a camera actually does when it looks at a flat surface from an angle, and which does not preserve parallelism. That's not a bug — it's exactly the phenomenon of a vanishing point, and it's the tool behind perspective correction and image stitching.

In [1]:

import numpy as np
import cv2
import matplotlib.pyplot as plt

Homogeneous coordinates¶

Represent a 2D point $(x, y)$ as a 3-vector $(x, y, 1)$. Any $3\times3$ matrix $H$ can then act on it by ordinary matrix multiplication; converting back to 2D means dividing by the third coordinate:

$$\begin{bmatrix}x'\\y'\\w'\end{bmatrix} = H\begin{bmatrix}x\\y\\1\end{bmatrix}, \qquad (x_{\text{2D}}', y_{\text{2D}}') = \left(\frac{x'}{w'}, \frac{y'}{w'}\right)$$

Two big payoffs: first, a $3\times3$ matrix can represent translation too (impossible with a $2\times2$ matrix alone, which is why Lesson 8 needed a separate translation vector $t$). Second, that division by $w'$ is what lets a homography model genuine perspective effects — not just the affine transforms of Lesson 8, but photos where parallel lines converge toward a vanishing point, exactly what a real camera does when it looks at a flat surface from an angle. Also note that $H$ and $cH$ (any nonzero scalar multiple) represent the exact same transform, since the division cancels the scale — a homography has only 8 independent degrees of freedom, not 9.

Parallel lines stop being parallel¶

The updated hierarchy, extending Lesson 8's table:

Transform	Matrix form	Preserves
Affine	$\begin{bmatrix}a&b&t_x\\c&d&t_y\\0&0&1\end{bmatrix}$	parallelism
Projective	$\begin{bmatrix}a&b&t_x\\c&d&t_y\\g&h&1\end{bmatrix}$	straight lines only

The only difference is a nonzero bottom row $(g, h)$ — and that alone is enough to destroy parallelism. We demonstrate directly: take two parallel horizontal lines and apply a homography with a nonzero bottom row.

In [2]:

H = np.array([
    [1,     0.2,   0],
    [0.1,   1,     0],
    [0.02,  0.015, 1],
], dtype=np.float64)

def apply_homography(points, H):
    homogeneous = np.hstack([points, np.ones((len(points), 1))])
    transformed = (H @ homogeneous.T).T
    return transformed[:, :2] / transformed[:, 2:3]

line1 = np.array([[0, 0], [1, 0]], dtype=np.float64)   # y = 0
line2 = np.array([[0, 1], [1, 1]], dtype=np.float64)   # y = 1, parallel to line1

l1_transformed = apply_homography(line1, H)
l2_transformed = apply_homography(line2, H)

print('line 1, transformed:', l1_transformed)
print('line 2, transformed:', l2_transformed)

def extended(p1, p2, length):
    direction = (p2 - p1)
    return p1 - length * direction, p2 + length * direction

fig, axes = plt.subplots(1, 2, figsize=(9, 4))

for line, color in [(line1, 'tab:blue'), (line2, 'tab:orange')]:
    a, b = extended(line[0], line[1], length=3)
    axes[0].plot([a[0], b[0]], [a[1], b[1]], color=color)
    axes[0].plot(line[:, 0], line[:, 1], 'o', color=color)
axes[0].set_xlim(-2, 4)
axes[0].set_ylim(-2, 3)
axes[0].set_title('Before: two parallel lines')

for line, color in [(l1_transformed, 'tab:blue'), (l2_transformed, 'tab:orange')]:
    a, b = extended(line[0], line[1], length=60)
    axes[1].plot([a[0], b[0]], [a[1], b[1]], color=color)
    axes[1].plot(line[:, 0], line[:, 1], 'o', color=color)
axes[1].set_xlim(-5, 55)
axes[1].set_ylim(-2, 8)
axes[1].set_title('After: a homography with a nonzero bottom row')

plt.tight_layout()
plt.show()

line 1, transformed: [[0.         0.        ]
 [0.98039216 0.09803922]]
line 2, transformed: [[0.19704433 0.98522167]
 [1.15942029 1.06280193]]

No description has been provided for this image

Vanishing points, computed two ways¶

Where do these two lines actually meet? We can find it two ways: intersecting the transformed line segments directly, or — more elegantly — transforming the point at infinity in the original lines' shared direction $(1, 0)$, represented in homogeneous coordinates as $(1, 0, 0)$ (a nonzero third coordinate would make it a finite point; zero means "infinitely far away"). Applying $H$ to that point at infinity should land exactly on the vanishing point.

In [3]:

def line_intersection(p1, p2, p3, p4):
    A = np.array([[p2[0] - p1[0], -(p4[0] - p3[0])],
                  [p2[1] - p1[1], -(p4[1] - p3[1])]])
    b = np.array([p3[0] - p1[0], p3[1] - p1[1]])
    t = np.linalg.solve(A, b)[0]
    return p1 + t * (p2 - p1)

vanishing_point_direct = line_intersection(l1_transformed[0], l1_transformed[1],
                                            l2_transformed[0], l2_transformed[1])

point_at_infinity = np.array([1.0, 0.0, 0.0])  # direction (1,0), infinitely far away
transformed_infinity = H @ point_at_infinity
vanishing_point_via_infinity = transformed_infinity[:2] / transformed_infinity[2]

print('vanishing point (line intersection): ', vanishing_point_direct)
print('vanishing point (H @ infinity trick):', vanishing_point_via_infinity)

vanishing point (line intersection):  [50.  5.]
vanishing point (H @ infinity trick): [50.  5.]

Lines in homogeneous coordinates¶

Points aren't the only thing homogeneous coordinates represent elegantly — a line $ax+by+c=0$ is just as naturally the 3-vector $l=(a,b,c)$, with a point $p=(x,y,1)$ lying on it exactly when the inner product $l\cdot p= l^\top p = 0$. Two useful consequences follow from the cross product: the line through two points is $l = p_1\times p_2$ (it's orthogonal to both, so both dot products vanish), and by the same logic in reverse, the intersection of two lines is $p = l_1\times l_2$. line_intersection above solved a $2\times2$ linear system by hand; the cross product does the same job in one line of code.

In [4]:

def line_through(p1, p2):
    return np.cross([p1[0], p1[1], 1.0], [p2[0], p2[1], 1.0])

def intersect(l1, l2):
    p = np.cross(l1, l2)
    return p[:2] / p[2]

line1_h = line_through(l1_transformed[0], l1_transformed[1])
line2_h = line_through(l2_transformed[0], l2_transformed[1])
vanishing_point_cross = intersect(line1_h, line2_h)

print('vanishing point (cross-product trick):', vanishing_point_cross)

vanishing point (cross-product trick): [50.  5.]

Matches the two earlier methods exactly. Points and lines (in the 2D plane) turn out to play completely symmetric roles in homogeneous coordinates — the same cross-product operation builds a line from two points or a point from two lines.

When does a homography actually apply?¶

A homography exactly relates two photos in exactly two situations: photographing a flat surface from two different positions (the case used below), or photographing any 3D scene, flat or not, from the same camera position while only rotating the camera between shots. That second case is what makes Lesson 24's panorama stitching work on ordinary, non-planar scenes — a translating camera or a genuinely 3D scene viewed from two different positions has no single homography that relates the two images exactly.

Perspective correction: rectifying a photographed document¶

The most common practical use of a homography: given 4 point correspondences (e.g. the 4 corners of a document, clicked by a user or found automatically), cv2.getPerspectiveTransform solves for the unique homography mapping one set of 4 points to the other, and cv2.warpPerspective applies it. The trapezoid shape a photographed rectangle takes on — edges that were parallel in real life visibly converging in the photo — is called keystoning, and it's exactly the vanishing-point effect above, now applied to a document instead of a pair of abstract lines. (Below, H_rectify is computed from the point correspondences in reverse order, but it's exactly np.linalg.inv(H_distort), since undoing a homography and inverting its matrix are the same operation.)

In [5]:

document = np.full((300, 220, 3), 255, dtype=np.uint8)
cv2.rectangle(document, (20, 20), (200, 280), (0, 0, 0), 3)
cv2.putText(document, 'HELLO', (30, 150), cv2.FONT_HERSHEY_SIMPLEX, 1.2, (0, 0, 0), 3)

h, w = document.shape[:2]
corners = np.float32([[0, 0], [w, 0], [w, h], [0, h]])
photographed_corners = np.float32([[45, 15], [175, 25], [210, 290], [5, 280]])  # strong keystone: top edge narrower than bottom

H_distort = cv2.getPerspectiveTransform(corners, photographed_corners)
photographed = cv2.warpPerspective(document, H_distort, (w, h))

H_rectify = cv2.getPerspectiveTransform(photographed_corners, corners)
rectified = cv2.warpPerspective(photographed, H_rectify, (w, h))

fig, axes = plt.subplots(1, 3, figsize=(9, 4))
for ax, im, title in zip(axes, [document, photographed, rectified],
                          ['Original document', 'Photographed at an angle\n(simulated)', 'Rectified\n(from 4 corner clicks)']):
    ax.imshow(cv2.cvtColor(im, cv2.COLOR_BGR2RGB))
    ax.set_title(title, fontsize=9)
    ax.axis('off')
plt.tight_layout()
plt.show()

error = np.abs(rectified.astype(int) - document.astype(int))
print(f'mean abs pixel difference after distort-then-rectify round trip: {error.mean():.2f}  (small resampling loss only)')

mean abs pixel difference after distort-then-rectify round trip: 4.59  (small resampling loss only)

Solving for a homography from many correspondences¶

4 correspondences exactly determine a homography's 8 degrees of freedom, with no slack for error. With more than 4 (typically from automatic feature matching, Lesson 19), the problem becomes an overdetermined least-squares fit — the Direct Linear Transform (DLT) algorithm that cv2.findHomography implements, optionally wrapped in RANSAC (Lesson 22) to reject bad correspondences (mismatched features) as outliers. cv2.findHomography also normalizes the point coordinates internally before solving, for numerical conditioning — the same trick made explicit in Lesson 26's 8-point algorithm.

In [6]:

rng = np.random.default_rng(0)
H_true = np.array([[1, 0.2, 10], [0.05, 1, 5], [0.0008, 0.0003, 1]])

pts1 = rng.uniform(0, 200, (30, 2))
pts2 = apply_homography(pts1, H_true)

H_estimated, inlier_mask = cv2.findHomography(pts1, pts2, cv2.RANSAC, 3.0)

print('true homography (scale-normalized):\n', np.round(H_true / H_true[2, 2], 4))
print('estimated homography:\n', np.round(H_estimated, 4))
print(f'inliers: {int(inlier_mask.sum())} / {len(inlier_mask)}')

true homography (scale-normalized):
 [[1.e+00 2.e-01 1.e+01]
 [5.e-02 1.e+00 5.e+00]
 [8.e-04 3.e-04 1.e+00]]
estimated homography:
 [[1.e+00 2.e-01 1.e+01]
 [5.e-02 1.e+00 5.e+00]
 [8.e-04 3.e-04 1.e+00]]
inliers: 30 / 30

With clean correspondences, findHomography recovers H_true almost exactly — this is exactly the last step of a typical image-stitching pipeline: detect and match SIFT features between two overlapping photos (Lesson 19), then solve for the homography that aligns one onto the other.

RANSAC in practice: homography estimation with bad matches¶

Even though the code above passes cv2.RANSAC to findHomography, with clean correspondences there is nothing for it to reject — it makes no visible difference. What effect does RANSAC actually have? Since real feature matching (Lesson 19) inevitably produces mismatched features, even after a ratio test, we need a way to discard such outliers. To illustrate this, we inject 15 deliberately garbage correspondences alongside the 30 genuine ones and compare the RANSAC fit against a plain least-squares fit on the exact same contaminated data.

In [7]:

pts1_bad = rng.uniform(0, 200, (15, 2))
pts2_bad = rng.uniform(0, 300, (15, 2))  # unrelated -- simulated bad matches

pts1_contaminated = np.vstack([pts1, pts1_bad])
pts2_contaminated = np.vstack([pts2, pts2_bad])

H_ransac, mask = cv2.findHomography(pts1_contaminated, pts2_contaminated, cv2.RANSAC, 3.0)
H_plain, _ = cv2.findHomography(pts1_contaminated, pts2_contaminated, 0)  # method=0: plain least squares, no outlier rejection

def mean_reprojection_error(H, pts1, pts2):
    return np.linalg.norm(apply_homography(pts1, H) - pts2, axis=1).mean()

print(f'inliers found by RANSAC: {int(mask.sum())} / {len(mask)}  (30 correspondences were genuinely correct)')
print()
print(f'mean reprojection error on the TRUE inliers:')
print(f'  RANSAC fit: {mean_reprojection_error(H_ransac, pts1, pts2):.2e} pixels')
print(f'  plain fit:  {mean_reprojection_error(H_plain, pts1, pts2):.2e} pixels')

inliers found by RANSAC: 30 / 45  (30 correspondences were genuinely correct)

mean reprojection error on the TRUE inliers:
  RANSAC fit: 1.49e-06 pixels
  plain fit:  1.06e+04 pixels

RANSAC finds exactly the 30 genuine correspondences and reconstructs the homography to essentially machine precision. The plain least-squares fit, given the exact same data, is off by many orders of magnitude more error — effectively useless for anything requiring pixel-level accuracy — because 15 bad matches out of 45 (33%) was more than enough to noticeably corrupt an unweighted sum-of-squares fit.

Looking ahead: projective geometry in 3D¶

Everything in this lesson has been 2D-to-2D: a homography relates two planes (images or flat surfaces in the world). Lesson 25 extends the same projective machinery one dimension further, deriving the camera's projection matrix that maps a 3D scene onto a 2D image in the first place.

Exercise¶

Add zero-mean Gaussian noise (e.g. std 2 pixels) to pts2 before calling cv2.findHomography. How much does the estimated homography drift from H_true, and does increasing the number of correspondences (say, from 30 to 200) reduce that drift?
Increase the number of injected bad matches (pts1_bad/pts2_bad) from 15 until RANSAC starts to fail. Using the iteration-count formula from Lesson 22 (minimal sample size 4 for a homography), roughly what outlier fraction is that, and how many iterations would it theoretically require to stay 99% confident of success?
In the vanishing-point demo, change the two lines to be vertical instead of horizontal (e.g. x=0 and x=1), and predict, then verify, where their vanishing point ends up using the point-at-infinity trick with direction (0, 1, 0).

← Lesson 22: Model Fitting, Robust Estimation, and RANSAC

Lesson 24: Image Stitching and Mosaicking →