Homography

A homography is the most general invertible map of the plane that keeps straight lines straight. In computer vision it is the transformation between two images of the same plane, such as a document, a floor, a building façade, or a calibration target, and between two images taken by a camera that rotates without moving its center. Because it is a $3 \times 3$ matrix that can be estimated from as few as four point correspondences, the homography is the workhorse of panorama stitching, planar rectification, augmented reality on flat surfaces, and camera calibration.

Definition

A homography, also called a projective transformation, projectivity, or collineation, is an invertible linear map of the projective plane $\mathbb{P}^2$, written in homogeneous coordinates as

$$\mathbf{x}' \sim H\mathbf{x}, \qquad H \in \mathbb{R}^{3 \times 3}, \quad \det H \neq 0,$$

where $\sim$ means equality up to a non-zero scale factor. Since $H$ and $kH$ define the same map, a homography has $9 - 1 = 8$ degrees of freedom. In Cartesian coordinates the map is a ratio of linear functions:

$$x' = \frac{h_{11}x + h_{12}y + h_{13}}{h_{31}x + h_{32}y + h_{33}}, \qquad y' = \frac{h_{21}x + h_{22}y + h_{23}}{h_{31}x + h_{32}y + h_{33}}.$$

The bottom row is what separates a homography from an affine map: when $h_{31} = h_{32} = 0$, the denominator is constant and parallel lines stay parallel. Projective geometry places the homography at the top of the hierarchy of planar transformations, above affine, similarity, and Euclidean maps, and describes what it preserves: collinearity, incidence, and the cross-ratio.

Intuition

Photograph a poster from the side. The rectangle becomes a general quadrilateral, the far edge shrinks, and the parallel edges converge, but every straight line on the poster is still straight in the photo. Any such picture of the poster can be turned into any other, or into a frontal view, by one homography, and four corners are enough to pin it down: each corner fixes two of the eight degrees of freedom.

The same holds for a camera that only turns. Rotating about the center of projection changes which ray each pixel sees but not which rays exist, so the new image is a re-arrangement of the old one, independent of how far away anything is. That is why a panorama can be stitched from photos taken while turning in place.

Formal Definition

When two images are related by a homography

A plane seen from two viewpoints. Let the first camera be $K[I \mid \mathbf{0}]$ and the second $K'[R \mid \mathbf{t}]$, and let the scene plane satisfy $\mathbf{n}^\top\mathbf{X} + d = 0$ in the first camera’s coordinates. Images of points on the plane are related by the plane-induced homography (Hartley and Zisserman, 2004, §13.1)

$$H = K'\left(R - \frac{\mathbf{t}\,\mathbf{n}^\top}{d}\right)K^{-1}.$$

Sign conventions differ between texts: writing the plane as $\mathbf{n}^\top\mathbf{X} = d$ instead, with $d$ the distance to the plane, changes the minus to a plus. The code below checks the form given here. The translation appears only as $\mathbf{t}/d$, so a homography cannot tell a small, close plane from a large, distant one.

A camera rotating about its center. With $\mathbf{t} = \mathbf{0}$ the plane drops out, and every scene point, at any depth, maps by

$$H = K' R K^{-1},$$

which is $KRK^{-1}$ when the intrinsic camera parameters do not change (Hartley and Zisserman, 2004, §8.4).

A plane and its image. A world plane $Z = 0$ maps into a single image by $H = K[\mathbf{r}_1\;\mathbf{r}_2\;\mathbf{t}]$, the camera matrix with its third column removed. This is the relation that planar camera calibration exploits (Zhang, 2000).

Estimation by the direct linear transform

Each correspondence $\mathbf{x}_i \leftrightarrow \mathbf{x}'_i$ requires $\mathbf{x}'_i \times H\mathbf{x}_i = \mathbf{0}$, which gives two independent equations that are linear in the nine entries $\mathbf{h}$ of $H$:

$$\begin{bmatrix} \mathbf{0}^\top & -w'_i\mathbf{x}_i^\top & y'_i\mathbf{x}_i^\top \\ w'_i\mathbf{x}_i^\top & \mathbf{0}^\top & -x'_i\mathbf{x}_i^\top \end{bmatrix}\mathbf{h} = \mathbf{0},$$

with $\mathbf{x}'_i = (x'_i, y'_i, w'_i)^\top$. Stacking $n \ge 4$ correspondences gives a $2n \times 9$ system $A\mathbf{h} = \mathbf{0}$. With exactly four points in general position the solution is the one-dimensional null space of $A$; with more, noisy points, the direct linear transform (DLT) takes the unit vector minimizing $\lVert A\mathbf{h}\rVert$, the right singular vector of $A$ with the smallest singular value (Hartley and Zisserman, 2004, §4.1). This is the homogeneous form of least squares.

Normalization

The DLT is not invariant to the choice of image coordinates. In pixel coordinates, entries of $A$ range from about 1 to about $10^6$, the matrix is badly conditioned, and the solution is sensitive to noise. Hartley (1997) showed for the eight-point algorithm that a simple normalization removes most of this sensitivity, and Hartley and Zisserman (2004, §4.4) present the normalized DLT for homographies as essential rather than optional: translate each image’s points so that their centroid is at the origin and scale them so that their average distance from it is $\sqrt{2}$, estimate $\tilde{H}$ from the normalized points with similarity transforms $T$ and $T'$, and undo the normalization with $H = T'^{-1}\tilde{H}T$.

Algebraic and geometric error

The DLT minimizes an algebraic error, a quantity with no direct geometric meaning. Better estimates minimize a geometric error in pixels: the transfer error $\sum_i d(\mathbf{x}'_i, H\mathbf{x}_i)^2$ in the second image, the symmetric transfer error in both, or, when both images have noise, the reprojection error, which also estimates corrected point positions $\hat{\mathbf{x}}_i$ with $\hat{\mathbf{x}}'_i = H\hat{\mathbf{x}}_i$ and is the maximum likelihood estimate under Gaussian noise. The Sampson error is a first-order approximation to the reprojection error that avoids estimating the corrected points. These costs are minimized iteratively, for example by Levenberg–Marquardt, starting from the normalized DLT (Hartley and Zisserman, 2004, §4.2–4.5).

Robust estimation

Correspondences from feature matching contain outliers, which ruin least squares. The standard pipeline runs the four-point DLT as the minimal solver inside RANSAC, scores each hypothesis by the transfer error of the matches, re-estimates $H$ from all inliers, and refines it by minimizing a geometric error (Hartley and Zisserman, 2004, §4.7–4.8). OpenCV’s findHomography follows this pattern with RANSAC, least median of squares, or the PROSAC-based RHO, discarding samples with collinear points and refining the result with Levenberg–Marquardt.

Properties

• Group structure. Homographies compose by matrix multiplication and invert by matrix inversion, so a chain of image-to-image homographies, as in a panorama, can be composed into one.
• Lines and conics. Points map by $H$, lines by $H^{-\top}$, and conics by $C' \sim H^{-\top}CH^{-1}$.
• Degeneracy. Four correspondences determine $H$ only if no three of the points are collinear in either image. If they are, $A$ has a null space of dimension greater than one and the estimate is undetermined or unstable; points spread over the image give a better-conditioned estimate than points clustered in one corner.
• Orientation. A homography estimated from points can map some points through the line at infinity, flipping them behind the camera. This cannot happen between two real images of points in front of the cameras, and it shows up as third coordinates of $H\mathbf{x}_i$ with different signs.

Decomposition into motion and structure

Given $K$, a plane-induced homography can be decomposed into the rotation $R$, the scaled translation $\mathbf{t}/d$, and the plane normal $\mathbf{n}$. Classical procedures compute the decomposition numerically; Malis and Vargas (2007) derived analytical expressions for the solutions as functions of the entries of $H$, which also expose how the solutions relate to each other. The decomposition is ambiguous: in general there are four solutions, two pairs related by a sign change, and requiring the observed points to lie in front of both cameras reduces them to two. Choosing between the last two needs further information, such as an approximate normal or a second plane. The scale of $\mathbf{t}$ is never recovered, and when the motion is a pure rotation, the plane normal is undefined. OpenCV implements the decomposition as decomposeHomographyMat and the visibility test as filterHomographyDecompByVisibleRefpoints.

Examples

The code below builds a plane-induced homography from a camera pair, checks it, and checks the pure-rotation formula. It then estimates $H$ with the DLT from ten noisy correspondences, with and without Hartley normalization, in a setting where all points lie within a patch about 160 pixels wide, and scores each estimate by its RMS transfer error on noise-free points of the patch.

import numpy as np
from scipy.spatial.transform import Rotation

rng = np.random.default_rng(1)


def to_h(p):
    return np.c_[p, np.ones(len(p))]


def apply(H, p):
    q = to_h(p) @ H.T
    return q[:, :2] / q[:, 2:]


# A camera K and a second camera [R | t]; the first camera is K [I | 0].
K = np.array([[800.0, 0, 640], [0, 800, 480], [0, 0, 1]])
R = Rotation.from_euler("xyz", [4, -12, 3], degrees=True).as_matrix()
t = np.array([0.6, -0.1, 0.2])

# Scene plane n^T X + d = 0 in the first camera's frame (Hartley & Zisserman's convention).
n = np.array([0.1, -0.2, -1.0]); n /= np.linalg.norm(n)
d = 5.0
H_true = K @ (R - np.outer(t, n) / d) @ np.linalg.inv(K)

# Check: points on the plane, projected by both cameras, are related by H_true.
uv = rng.uniform(-0.5, 0.5, (200, 2))      # a patch of the plane, about 160 px wide in view 1
X = np.c_[uv, np.zeros(200)]
X[:, 2] = (-d - X[:, :2] @ n[:2]) / n[2]          # solve n^T X + d = 0 for Z
x1 = apply(K, X[:, :2] / X[:, 2:])                  # K [I | 0] X
Xc2 = X @ R.T + t
x2 = (Xc2 @ K.T)[:, :2] / (Xc2 @ K.T)[:, 2:]
print("plane-induced H maps view 1 to view 2, max error (px): %.1e"
      % np.abs(apply(H_true, x1) - x2).max())

# Pure rotation: H = K R K^-1 for any scene depth.
Xr = rng.uniform(-3, 3, (50, 3)); Xr[:, 2] = rng.uniform(2, 50, 50)
y1 = (Xr @ K.T)[:, :2] / (Xr @ K.T)[:, 2:]
y2 = ((Xr @ R.T) @ K.T)[:, :2] / ((Xr @ R.T) @ K.T)[:, 2:]
print("rotation-only H = K R K^-1, max error (px): %.1e"
      % np.abs(apply(K @ R @ np.linalg.inv(K), y1) - y2).max())


def normalize(p):
    """Similarity T: centroid to the origin, mean distance from it sqrt(2)."""
    c = p.mean(0)
    s = np.sqrt(2) / np.linalg.norm(p - c, axis=1).mean()
    return np.array([[s, 0, -s * c[0]], [0, s, -s * c[1]], [0, 0, 1]])


def dlt(p, q, normalized=True):
    """Estimate H with q ~ H p from N >= 4 correspondences."""
    T1 = normalize(p) if normalized else np.eye(3)
    T2 = normalize(q) if normalized else np.eye(3)
    a, b = apply(T1, p), apply(T2, q)
    rows = []
    for (x, y), (u, v) in zip(a, b):
        rows.append([0, 0, 0, -x, -y, -1, v * x, v * y, v])
        rows.append([x, y, 1, 0, 0, 0, -u * x, -u * y, -u])
    _, _, Vt = np.linalg.svd(np.array(rows))
    Hn = Vt[-1].reshape(3, 3)                        # null vector of A, ||h|| = 1
    H = np.linalg.inv(T2) @ Hn @ T1                  # undo the normalization
    return H / H[2, 2]


# Exact recovery from the minimal 4 points without noise.
H4 = dlt(x1[:4], x2[:4])
print("4 exact points: max |H - H_true| = %.1e" % np.abs(H4 - H_true / H_true[2, 2]).max())


def design_cond(p, q):
    """Condition number (largest / second-smallest singular value) of the DLT matrix."""
    rows = []
    for (x, y), (u, v) in zip(p, q):
        rows += [[0, 0, 0, -x, -y, -1, v * x, v * y, v], [x, y, 1, 0, 0, 0, -u * x, -u * y, -u]]
    s = np.linalg.svd(np.array(rows), compute_uv=False)
    return s[0] / s[-2]


print("DLT matrix condition: pixels %.1e, normalized %.1e"
      % (design_cond(x1[:10], x2[:10]),
         design_cond(apply(normalize(x1[:10]), x1[:10]), apply(normalize(x2[:10]), x2[:10]))))

# Noisy correspondences: 10 points, sigma = 1 px, 1000 trials. Score each H by its
# RMS transfer error on all 200 noise-free points of the patch.
sigma, N, trials = 1.0, 10, 1000
err = {True: [], False: []}
for _ in range(trials):
    idx = rng.choice(len(x1), N, replace=False)
    p = x1[idx] + sigma * rng.standard_normal((N, 2))
    q = x2[idx] + sigma * rng.standard_normal((N, 2))
    for norm in (True, False):
        H = dlt(p, q, norm)
        err[norm].append(np.sqrt(np.mean(np.sum((apply(H, x1) - x2) ** 2, 1))))
for norm, name in [(True, "normalized  "), (False, "unnormalized")]:
    e = np.array(err[norm])
    print(f"{name} DLT: RMS error median {np.median(e):6.2f} px, "
          f"90th percentile {np.percentile(e, 90):7.2f} px")

Output:

plane-induced H maps view 1 to view 2, max error (px): 2.3e-13
rotation-only H = K R K^-1, max error (px): 2.3e-13
4 exact points: max |H - H_true| = 1.0e-11
DLT matrix condition: pixels 6.6e+07, normalized 3.7e+00
normalized   DLT: RMS error median   1.55 px, 90th percentile    2.59 px
unnormalized DLT: RMS error median   6.11 px, 90th percentile  124.64 px

Both formulas reproduce the projections to rounding error, and four exact correspondences recover $H$. In pixel coordinates the DLT matrix has a condition number of nearly $10^8$; after normalization it is below 4. With one-pixel noise, the normalized DLT’s median error over 1,000 trials is 1.55 pixels, while the unnormalized DLT’s is about four times larger and its worst tenth of estimates are off by more than 100 pixels. The effect depends on the configuration: with points spread across the whole image, the gap is much smaller, but normalization costs almost nothing.

Common Misconceptions

• “Any two views of a scene are related by a homography.” Only if the scene is planar or the camera does not translate. For a general 3D scene and a moving camera, points at different depths shift by different amounts (parallax), and the two views are related by epipolar geometry, which constrains a point only to a line. A homography fitted to such data aligns one plane, or the distant background, and misaligns everything else.
• “Four matches are enough in practice.” Four is the minimum, used inside RANSAC. The final estimate should use all inliers and a geometric error.
• “The decomposition recovers the camera motion.” It recovers it up to a two-fold ambiguity and an unknown scale, and only with known intrinsics.

Where It Is Used

• Image stitching and panoramas. Images taken by a rotating camera, or of a distant scene, are aligned by homographies estimated from matched features and then warped onto a common surface and blended (Szeliski, 2007). For a purely rotating camera with known $K$, parameterizing the alignment by rotations rather than general homographies reduces the unknowns.
• Planar rectification. A homography from four known points, such as the corners of a document, a sign, or a sports field, produces a fronto-parallel view in which lengths and angles on the plane can be measured.
• Camera calibration. Planar camera calibration estimates one homography per view of a flat target; each gives two constraints on the intrinsics (Zhang, 2000).
• Augmented reality and planar tracking. Tracking a planar marker or textured surface frame by frame amounts to estimating a homography, from which the pose of the plane follows by decomposition when $K$ is known.
• Geometric verification. In image retrieval and matching, a homography fitted with RANSAC checks whether candidate matches between two images are consistent.
• Learned matching. Synthetic homographic warps of real images supply training pairs with known correspondences, as in the homographic adaptation used to train SuperPoint, and homography estimation accuracy is a common way to evaluate matchers.