Direct Linear Transform (DLT)

The direct linear transform (DLT) estimates a projective mapping from point correspondences without iteration. Each correspondence gives equations that are linear in the unknown matrix entries; stacking them yields a homogeneous system $A\mathbf{m} = \mathbf{0}$ , and the singular value decomposition of $A$ gives the answer in one step. The same recipe estimates a homography from matched image points, a camera matrix from known 3D points and their images, and, in slightly different forms, a triangulated point or a fundamental matrix. In practice the DLT rarely gives the final answer. It is the fast, closed-form starting point that a geometric refinement then improves, and the minimal solver run inside RANSAC.

Problem

Many relations in multiple-view geometry have the form $\mathbf{x}_i \sim M\mathbf{X}_i$ : a known point $\mathbf{X}_i$ , in homogeneous coordinates, maps by an unknown matrix $M$ to a measured image point $\mathbf{x}_i$ , up to an unknown scale factor that differs for each point. For a homography, $\mathbf{X}_i$ is a point in another image and $M$ is $3 \times 3$ . For camera resectioning, $\mathbf{X}_i$ is a 3D point and $M$ is the $3 \times 4$ camera matrix $P$ of the pinhole camera model. Dividing out the per-point scale produces ratios, nonlinear in $M$ ; the DLT removes the scale without division.

Inputs and Outputs

Input: $n$ correspondences $\mathbf{X}_i \leftrightarrow \mathbf{x}_i$ (at least 4 for a homography, 6 for a camera matrix), with $\mathbf{X}_i$ 2D or 3D and $\mathbf{x}_i$ in pixels.
Output: the matrix $M$ , up to scale, normalized to unit Frobenius norm; optionally the singular values of $A$ , which show whether the data determined it.

Intuition

Two homogeneous vectors represent the same point exactly when they are parallel, that is, when their cross product is zero. Writing “ $\mathbf{x}_i$ is parallel to $M\mathbf{X}_i$ ” this way gives equations linear in the entries of $M$ , because $\mathbf{x}_i$ is known. Since $M$ and $2M$ describe the same mapping, the solution is a direction: the one that comes closest to satisfying all the equations.

Algorithm

Normalize each point set: translate it so its centroid is at the origin and scale it so the average distance from the origin is $\sqrt{2}$ for 2D points ( $\sqrt{3}$ for 3D points in the code below; some texts use the root-mean-square distance instead), giving similarity transforms $T$ (for $\mathbf{X}$ ) and $U$ (for $\mathbf{x}$ ).
Build $A$ : for each normalized correspondence, append the two rows derived below.
Solve: compute the SVD $A = W\Sigma V^\top$ and take $\tilde{\mathbf{m}}$ , the right singular vector of the smallest singular value; reshape it into $\tilde{M}$ .
Denormalize: $M = U^{-1}\tilde{M}\,T$ .
Refine (recommended): minimize a geometric error, such as the reprojection error, starting from $M$ .

Mathematical Formulation

Two equations per point

Let $\mathbf{x}_i = (x_i, y_i, w_i)^\top$ and write $\mathbf{m}^{j\top}$ for the $j$ -th row of $M$ , so that $M\mathbf{X}_i = (\mathbf{m}^{1\top}\mathbf{X}_i,\, \mathbf{m}^{2\top}\mathbf{X}_i,\, \mathbf{m}^{3\top}\mathbf{X}_i)^\top$ . The condition $\mathbf{x}_i \times M\mathbf{X}_i = \mathbf{0}$ is three equations, linear in the stacked rows $\mathbf{m} = (\mathbf{m}^1; \mathbf{m}^2; \mathbf{m}^3)$ :

\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 \\ -y_i\mathbf{X}_i^\top & x_i\mathbf{X}_i^\top & \mathbf{0}^\top \end{bmatrix}\mathbf{m} = \mathbf{0}.

Only two of them are independent. Multiplying the first row by $x_i$ and the second by $y_i$ and adding gives $-w_i$ times the third, so for $w_i \neq 0$ the third row adds nothing and is dropped. Geometrically, an image point fixes two numbers, its two coordinates, so it can constrain only two degrees of freedom.

The same block serves both problems. For a homography, $\mathbf{X}_i$ is a 3-vector and $A$ is $2n \times 9$ , exactly the system the homography article writes out. For a camera matrix, $\mathbf{X}_i = (X_i, Y_i, Z_i, 1)^\top$ is a 4-vector and $A$ is $2n \times 12$ .

Counting points

A homography has $9 - 1 = 8$ degrees of freedom, so four points in general position give eight equations and determine it. A camera matrix has $12 - 1 = 11$ , so five points give only ten equations and five and a half are needed: six points, with one equation to spare. With exact data in general position, the null space of $A$ is one-dimensional and the solution is exact. With noisy data and more equations than degrees of freedom, as already for six points and a camera, $A$ has full column rank and no exact solution exists.

The SVD solution and the norm constraint

Since $\mathbf{m}$ is defined only up to scale, $\mathbf{m} = \mathbf{0}$ must be excluded, and the DLT fixes the scale with $\lVert \mathbf{m} \rVert = 1$ :

\hat{\mathbf{m}} = \arg\min_{\lVert \mathbf{m} \rVert = 1} \lVert A\mathbf{m} \rVert,

whose solution is the right singular vector of $A$ for its smallest singular value. The singular value decomposition article proves this and explains why the gap between the two smallest singular values governs how far noise can move the answer. This is the homogeneous form of least squares.

The original photogrammetric formulation fixed the scale differently. Abdel-Aziz and Karara (1971) wrote the image coordinates as ratios of linear functions with eleven coefficients $L_1, \dots, L_{11}$ and the constant 1 in the denominator (in the form given by El-Ashmawy, 2018), which amounts to setting the last entry of $P$ to 1 and solving an ordinary least-squares problem. It follows from $P(0, 0, 0, 1)^\top = \mathbf{p}_4$ , the fourth column of $P$ , that this entry is proportional to the depth of the world origin, so the inhomogeneous form fails when the world origin lies near the plane through the camera center parallel to the image (our reasoning). The unit-norm constraint has no such blind spot.

Normalization

Rows of $A$ contain products such as $x_i X_i$ . With pixel coordinates in the hundreds and world coordinates in the thousands of millimetres, entries of one row range from 1 to about $10^6$ and the matrix is badly conditioned. Hartley (1997) showed for the eight-point algorithm that centering and scaling each point set before solving greatly improves accuracy, and attributed the improvement to better conditioning of the matrix in the eigenvalue problem (as summarized by Chojnacki et al., 2003; see the article on Hartley’s paper). Chojnacki et al. (2003) gave a second explanation: the normalized algorithm minimizes a different cost function whose terms are more balanced in spread, and therefore statistically better founded, than those of the unnormalized one.

A simpler reason (our reasoning): the constraint $\lVert \mathbf{m} \rVert = 1$ is not preserved by a change of coordinates, so without normalization the answer depends on the origin and units of both point sets, such as pixels counted from the image corner or the center, or millimetres versus metres. Normalizing both sets by a fixed rule removes that dependence. Afterwards, since $U\mathbf{x} \sim \tilde{M}\,T\mathbf{X}$ , the estimate in original coordinates is $M = U^{-1}\tilde{M}\,T$ .

Algebraic and geometric error

The DLT minimizes $\lVert A\mathbf{m} \rVert$ , an algebraic error: a sum of residuals that weight points unequally and have no unit of pixels. What should be minimized when the image measurements are noisy is a geometric error, the distance in pixels between measured and predicted points, whose minimum under Gaussian noise is the maximum-likelihood estimate (see reprojection error). That error is nonlinear in $\mathbf{m}$ and is minimized iteratively, for example by Levenberg–Marquardt as described in nonlinear least squares, with the normalized DLT as the starting point.

Degenerate Configurations

Some point sets do not determine $M$ however many points they contain. Both classic cases have the same structure (our derivation, confirmed by the code below):

Collinear points, homography. If all $\mathbf{X}_i$ lie on a line $\mathbf{l}$ , so that $\mathbf{l}^\top\mathbf{X}_i = 0$ , then $H + \mathbf{v}\mathbf{l}^\top$ maps every point exactly as $H$ does, for any 3-vector $\mathbf{v}$ . The null space of $A$ grows from one dimension to four. A minimal sample of four points with three on a line is degenerate too: its $8 \times 9$ matrix has a two-dimensional null space.
Coplanar points, camera matrix. If all $\mathbf{X}_i$ lie on a plane $\boldsymbol{\pi}$ , with $\boldsymbol{\pi}^\top\mathbf{X}_i = 0$ , then $P + \mathbf{v}\boldsymbol{\pi}^\top$ projects every point exactly as $P$ does, and again the null space is four-dimensional. Points on one plane determine only the homography from that plane to the image, eight of the eleven degrees of freedom. This is why classical calibration with the DLT uses a 3D object, while planar camera calibration estimates one homography per view and combines several views.

Nearly degenerate data are also harmful: points almost on a line, or a 3D object that is nearly flat, leave a second singular value close to the smallest, and noise then mixes their singular vectors.

Complexity

Building $A$ and computing the SVD of a $2n \times k$ matrix, with $k$ = 9 or 12, both take time linear in $n$ , so a minimal DLT can run thousands of times inside RANSAC.

Use Inside RANSAC

Correspondences from feature matching contain wrong matches that ruin a least-squares fit. RANSAC repeatedly fits $M$ to a random minimal sample, four matches for a homography or six 3D–2D pairs for a camera matrix, and keeps the hypothesis that most correspondences agree with. The DLT is a natural minimal solver: it is closed-form and cheap, and its singular values flag degenerate samples, such as collinear points, before they are scored. The final model is the normalized DLT on all inliers followed by geometric refinement. The perspective-n-point article runs this pipeline, with the calibrated DLT, on simulated data with about 40% outliers.

Implementation

The NumPy code below implements one DLT for both cases: dlt returns a homography for 2D points and a camera matrix for 3D points. The 54 inner corners of a checkerboard are seen in two $1280 \times 960$ views related by a known homography. A calibration object made of two orthogonal 300 mm faces, each carrying a $6 \times 6$ grid of targets, has its corner at (1.5, 0.8, 2.0) m in a room coordinate frame and is imaged by a camera with a 1,000-pixel focal length, once from about 0.7 m and once from about 3 m. Image points receive 1 pixel of Gaussian noise, over ten seeds, and each estimate is scored by the RMS distance between where it and the true model map the noise-free points. The script then counts the near-zero singular values of $A$ for two degenerate configurations.

import numpy as np


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


def project(M, p):
    """Apply a projective map (homography, camera, or normalizing T) to points."""
    q = hom(p) @ M.T
    return q[:, :-1] / q[:, -1:]


def normalizer(p):
    """Hartley normalization: centroid to the origin, mean distance sqrt(d)."""
    d = p.shape[1]
    c = p.mean(0)
    s = np.sqrt(d) / np.linalg.norm(p - c, axis=1).mean()
    T = np.eye(d + 1)
    T[:d, :d] *= s
    T[:d, d] = -s * c
    return T


def dlt_matrix(X, x):
    """Two rows per correspondence x ~ M X, from the cross product x × (M X) = 0."""
    rows = []
    for Xi, (u, v) in zip(hom(X), x):
        z = np.zeros_like(Xi)
        rows.append(np.r_[z, -Xi, v * Xi])
        rows.append(np.r_[Xi, z, -u * Xi])
    return np.array(rows)


def dlt(X, x, normalize=True):
    """Estimate M with x ~ M X: a 3x3 homography if X is 2D, a 3x4 camera if 3D."""
    T = normalizer(X) if normalize else np.eye(X.shape[1] + 1)
    U = normalizer(x) if normalize else np.eye(3)
    A = dlt_matrix(project(T, X), project(U, x))
    _, s, Vt = np.linalg.svd(A)
    M = Vt[-1].reshape(3, -1)                  # unit vector minimizing ||A m||
    M = np.linalg.inv(U) @ M @ T               # undo the normalization
    return M / np.linalg.norm(M), s[0] / s[-2]  # cond: sigma_1 / sigma_(n-1)


def rms(M, X, x_true):
    return np.sqrt(np.mean(np.sum((project(M, X) - x_true) ** 2, 1)))


# Homography: the 54 inner corners of a checkerboard seen in two 1280 x 960 views.
g = np.mgrid[0:9, 0:6].T.reshape(-1, 2) * 84.0 + [180, 150]   # view 1 (pixels)
H = np.array([[1.05, 0.12, 40.0], [-0.08, 0.95, 30.0], [1.2e-4, -2.2e-4, 1.0]])
g2 = project(H, g)                                            # view 2 (pixels)

# Camera: f = 1000 px, 1280 x 960 image. A calibration object with two orthogonal
# 300 mm faces, a 6 x 6 grid of targets on each, has its corner at (1.5, 0.8, 2) m.
K = np.array([[1000.0, 0, 640], [0, 1000, 480], [0, 0, 1]])
a = np.deg2rad(35)
R = np.array([[np.cos(a), 0, -np.sin(a)], [0, 1, 0], [np.sin(a), 0, np.cos(a)]])
u = np.mgrid[0:6, 0:6].T.reshape(-1, 2) * 60.0
Xo = np.r_[np.c_[u, np.zeros(36)], np.c_[np.zeros(36), u[:, 1], u[:, 0]]]
origin = np.array([1500.0, 800, 2000])                 # object corner in room frame (mm)
X = Xo + origin
P_near = K @ np.c_[R, R @ -origin + [-60, -150, 700]]  # object about 0.7 m away
P_far = K @ np.c_[R, R @ -origin + [-60, -150, 3000]]  # object about 3 m away
x, x_far = project(P_near, X), project(P_far, X)
for name, p in [("board", g2), ("object, near", x), ("object, far", x_far)]:
    w, h = p.max(0) - p.min(0)
    print(f"{name}: spans {w:.0f} x {h:.0f} px")

# 1. Minimal, exact data: 4 points fix H, 6 non-coplanar points fix P.
h4, p6 = [0, 8, 45, 53], [0, 5, 35, 41, 50, 71]
print("exact minimal sets: H from 4 points %.1e px, P from 6 points %.1e px"
      % (rms(dlt(g[h4], g2[h4])[0], g, g2), rms(dlt(X[p6], x[p6])[0], X, x)))

# 2. Noise of 1 px, ten seeds. Error: RMS distance between the projections of the
# noise-free points by the estimate and by the true model.
h8, p8 = [0, 4, 8, 22, 31, 45, 49, 53], [0, 5, 14, 30, 35, 41, 50, 71]
cases = [("H, 54 pts       ", g, g2, slice(None)), ("H, 8 pts        ", g, g2, h8),
         ("P near, 72 pts  ", X, x, slice(None)), ("P near, 8 pts   ", X, x, p8),
         ("P far, 72 pts   ", X, x_far, slice(None)), ("P far, 8 pts    ", X, x_far, p8)]
for name, A3, a2, idx in cases:
    for norm in (False, True):
        c, e = [], []
        for seed in range(10):
            noisy = a2 + np.random.default_rng(seed).normal(0, 1.0, a2.shape)
            M, k = dlt(A3[idx], noisy[idx], norm)
            c.append(k)
            e.append(rms(M, A3, a2))
        print(f"{name} {'normalized  ' if norm else 'unnormalized'}: cond(A) "
              f"{min(c):.1e}-{max(c):.1e}  error median {np.median(e):5.2f} px, "
              f"max {max(e):6.2f} px")


# 3. Degenerate configurations, exact data: count singular values near zero.
def null_dim(A3, a2):
    A = dlt_matrix(project(normalizer(A3), A3), project(normalizer(a2), a2))
    s = np.linalg.svd(A, compute_uv=False)
    return A.shape[1] - np.sum(s > 1e-9 * s[0]), s[-5:] / s[0]


for name, A3, a2, idx in [("H, 4 corners, 3 on a row ", g, g2, np.array([0, 1, 2, 53])),
                          ("H, 9 corners on one row  ", g, g2, np.arange(9)),
                          ("H, all 54 corners        ", g, g2, np.arange(54)),
                          ("P, 36 targets on one face", X, x, np.arange(36)),
                          ("P, all 72 targets        ", X, x, np.arange(72))]:
    k, s = null_dim(A3[idx], a2[idx])
    print(f"{name}: null space dim {k}, smallest s / s1 {np.array2string(s, precision=1)}")

Output:

board: spans 777 x 524 px
object, near: spans 458 x 429 px
object, far: spans 130 x 100 px
exact minimal sets: H from 4 points 1.7e-13 px, P from 6 points 6.0e-13 px
H, 54 pts        unnormalized: cond(A) 2.2e+06-2.3e+06  error median  0.48 px, max   0.60 px
H, 54 pts        normalized  : cond(A) 3.4e+00-3.4e+00  error median  0.34 px, max   0.50 px
H, 8 pts         unnormalized: cond(A) 2.0e+06-2.1e+06  error median  0.79 px, max   1.34 px
H, 8 pts         normalized  : cond(A) 3.6e+00-3.6e+00  error median  0.75 px, max   1.33 px
P near, 72 pts   unnormalized: cond(A) 1.3e+08-1.4e+08  error median  0.39 px, max   0.54 px
P near, 72 pts   normalized  : cond(A) 5.8e+00-5.8e+00  error median  0.39 px, max   0.50 px
P near, 8 pts    unnormalized: cond(A) 9.6e+07-9.8e+07  error median  0.84 px, max   1.16 px
P near, 8 pts    normalized  : cond(A) 5.0e+00-5.0e+00  error median  0.84 px, max   1.16 px
P far, 72 pts    unnormalized: cond(A) 3.2e+08-4.0e+08  error median 86.02 px, max 449.88 px
P far, 72 pts    normalized  : cond(A) 6.0e+00-6.1e+00  error median  0.39 px, max   0.50 px
P far, 8 pts     unnormalized: cond(A) 3.2e+08-9.1e+08  error median  1.14 px, max   9.64 px
P far, 8 pts     normalized  : cond(A) 5.1e+00-5.2e+00  error median  0.81 px, max   1.11 px
H, 4 corners, 3 on a row : null space dim 2, smallest s / s1 [2.3e-01 8.3e-02 4.8e-02 7.3e-03 3.1e-17]
H, 9 corners on one row  : null space dim 4, smallest s / s1 [1.8e-01 5.8e-17 6.7e-19 3.3e-19 0.0e+00]
H, all 54 corners        : null space dim 1, smallest s / s1 [3.9e-01 3.4e-01 3.2e-01 2.9e-01 1.1e-16]
P, 36 targets on one face: null space dim 4, smallest s / s1 [2.9e-01 5.7e-16 5.6e-17 8.9e-18 0.0e+00]
P, all 72 targets        : null space dim 1, smallest s / s1 [3.1e-01 3.0e-01 2.9e-01 1.7e-01 4.9e-16]

With exact data, both minimal problems are solved to rounding error. Without normalization, the condition number of $A$ (the ratio of the largest to the second-smallest singular value, since the smallest is ideally zero) is about $10^6$ for the homography and $10^8$ to $10^9$ for the camera, where world coordinates in millimetres meet pixel coordinates; with normalization it is below 7 in every case.

Bad conditioning does not always mean a bad answer. When the points fill much of the image, the unnormalized estimates are somewhat worse for the homography (median 0.48 against 0.34 pixels with all points) and no worse for the camera seen from 0.7 m. The penalty appears when the object covers a small part of the image: seen from 3 m, it spans about 130 by 100 pixels, and the unnormalized estimate from all 72 points has a median error of 86 pixels, while the normalized one stays at 0.39 pixels, the same as from close up. Adding points made the unnormalized estimate worse in this case, not better, which suggests a systematic error of the unbalanced cost rather than noise that averages out (our reading). The homography article shows the same effect for points in a small image patch. Normalization costs two $3 \times 3$ or $4 \times 4$ matrix products, so there is no reason to skip it.

The degenerate cases behave as derived above. Nine corners on one row of the board, and the 36 targets of one face, each leave a four-dimensional null space, in which the last singular vector is an arbitrary member; four corners with three on a row leave a two-dimensional one. The full point sets leave exactly one.

Limitations

Algebraic, not geometric, error. The estimate is not the maximum-likelihood one and should be refined when accuracy matters; in exchange it needs no starting guess and has no local minima.
No constraints on structure. The camera DLT estimates a general $3 \times 4$ matrix: it cannot impose known intrinsics, zero skew, or a rotation that is exactly orthonormal. When the intrinsics are known, perspective-n-point solvers use them instead of estimating them.
No lens distortion. The model is projective. Radial distortion must be removed first or estimated in the nonlinear refinement.

Variants

Modified DLT. El-Ashmawy (2018) states that the eleven DLT coefficients contain only ten independent parameters, and describes the modified DLT, which adds the missing relation between them as a nonlinear constraint. In the projective camera model the eleventh degree of freedom is the skew of the image axes, so the constraint amounts to requiring orthogonal axes (our reasoning).
Linear triangulation. With known cameras and an unknown 3D point, the same cross-product construction gives a $2n \times 4$ system; see triangulation.
Eight-point algorithm. The epipolar constraint gives one equation per correspondence, linear in the nine entries of the fundamental matrix, solved by the same SVD step plus a rank-2 correction.
Linear pose. With known intrinsics, the DLT on calibrated coordinates followed by projection onto the nearest rotation gives a linear camera pose, the simplest non-minimal PnP method.

History

Abdel-Aziz and Karara introduced the direct linear transformation to photogrammetry at a 1971 symposium on close-range photogrammetry in Urbana, Illinois (El-Ashmawy, 2018). A paper by the same authors with nearly the same title appeared in Photogrammetric Engineering & Remote Sensing in 2015, evidently a reprint; El-Ashmawy (2018) cites both. The 2015 title speaks of comparator coordinates, measured on the photograph, and the method related them to object-space coordinates through the eleven coefficients directly, without requiring the interior orientation of the camera to be known in advance, which the resection methods based on the collinearity equations need (El-Ashmawy, 2018). El-Ashmawy (2018) lists what made it popular in close-range photogrammetry, computer vision, robotics and biomechanics: the image coordinates may be in a non-orthogonal system with unequal scales, the origin of the coordinate system is arbitrary, and the principal distance may be unknown and may change from image to image. In computer vision the method is usually written in the homogeneous form above and solved by the SVD after Hartley’s normalization.