Triangulation & 3D Point Recovery

"We were supposed to meet at the point. I came in from the left camera, she came in from the right, and we missed each other by half a millimeter, forever. Now a least-squares solver decides where our relationship 'really' was."

A Slightly Skew Triangulation Ray

Section 13.5 recovered depth densely, but only for the rectified, calibrated, short-baseline world of stereo rigs. The general problem is sparser and wilder: two arbitrary cameras, poses recovered by Section 13.2's essential-matrix machinery, and a few thousand verified matches from Chapter 10. Each match nominates two rays in space; geometry says they should cross at the scene point. Then measurement noise has its say.

1. Why the Rays Miss Beginner

Each detected keypoint is off by some fraction of a pixel (recall the sub-pixel refinement battles of Chapter 10); each camera pose carries its own estimation error. Back-project two such corrupted measurements and you get two lines in 3D space, and two lines in 3D, unlike two lines in a plane, generically do not intersect at all: they pass by each other, skew, at some closest-approach distance. "The triangulated point" is therefore a decision, not an observation, and Figure 13.6.1 shows both the decision and the geometry that controls its quality: the angle between the rays.

The right panel of Figure 13.6.1 is the section's recurring warning. When the parallax angle between rays is small (short baseline relative to depth), pixel noise barely changes where the rays are but enormously changes where they nearly cross: the uncertainty region stretches into a sliver along the viewing direction. This is Section 13.5's quadratic law wearing general-geometry clothes, and it is why every serious pipeline filters triangulated points by parallax angle before trusting them.

2. The Midpoint Method Intermediate

The estimate drawn in Figure 13.6.1's left panel is the midpoint method: find the unique points $\mathbf{p}_1 = \mathbf{C}_L + s\,\mathbf{r}_L$ and $\mathbf{p}_2 = \mathbf{C}_R + u\,\mathbf{r}_R$ minimizing $\lVert \mathbf{p}_1 - \mathbf{p}_2 \rVert$ (two linear equations from perpendicularity of the connecting segment to both rays), and return their midpoint. It is intuitive, fast, and fine for visualization or initialization. Its weakness is principle: it minimizes a distance in 3D space, where errors are anisotropic and depth-dependent (that sliver again), so it silently weights distant points' wild depth axes the same as their well-measured lateral axes. The measurements we actually trust live in the image plane, which suggests minimizing there instead.

3. Linear Triangulation: The DLT Again Advanced

The workhorse algorithm is the linear one, and it recycles this chapter's favorite move (one homogeneous equation set, one SVD null space). A camera with $3 \times 4$ projection matrix $P$ maps the homogeneous point $\mathbf{X}$ to $\mathbf{x} \sim P\mathbf{X}$; equality up to scale means the cross product vanishes, $\mathbf{x} \times P\mathbf{X} = \mathbf{0}$, giving two independent linear equations per view. Stack the four equations from two views into $A\mathbf{X} = \mathbf{0}$ and take the SVD null vector, exactly as in the eight-point algorithm of Section 13.2 and the homography DLT of Section 13.3:

import cv2
import numpy as np

def triangulate_dlt(P1, P2, x1, x2):
    """One 3D point from one match. P1, P2: 3x4. x1, x2: pixel coords."""
    A = np.stack([
        x1[0] * P1[2] - P1[0],      # x * (row 3) - (row 1) = 0
        x1[1] * P1[2] - P1[1],      # y * (row 3) - (row 2) = 0
        x2[0] * P2[2] - P2[0],
        x2[1] * P2[2] - P2[1]])
    _, _, Vt = np.linalg.svd(A)
    X = Vt[-1]
    return X[:3] / X[3]             # dehomogenize

# left camera at origin, right camera from Section 13.2's recoverPose
P1 = K @ np.hstack([np.eye(3), np.zeros((3, 1))])
P2 = K @ np.hstack([R, t.reshape(3, 1)])

X = triangulate_dlt(P1, P2, ptsL[0], ptsR[0])
print(np.round(X, 3))               # e.g. [-0.214  0.078  3.642]

# the same, vectorized over all matches, via OpenCV
Xh = cv2.triangulatePoints(P1, P2, ptsL.T.astype(np.float64),
                                   ptsR.T.astype(np.float64))
pts3d = (Xh[:3] / Xh[3]).T          # Nx3

Caveats that bite in practice: cv2.triangulatePoints wants points as $2 \times N$ arrays and returns homogeneous $4 \times N$ output, and forgetting the divide-by-$w$ is a rite of passage. The linear solve minimizes an algebraic quantity, not reprojection error; with Hartley normalization of image coordinates it is typically within a whisker of optimal at sane parallax, and a few Gauss-Newton steps on the reprojection objective (or Hartley-Sturm's exact polynomial method for the two-view case) close the remaining gap when millimeters matter. Because the result is computed in homogeneous coordinates, a point can come back with $w$ of the sign that places it behind a camera, which is not a bug but a verdict: that match was wrong, and the check belongs in the pipeline.

4. The Complete Two-View Reconstructor Intermediate

Everything this chapter built now snaps together into a single pipeline: matches (Chapter 10), essential matrix and pose (Section 13.2), triangulation (this section), and the two filters this section has motivated: cheirality (positive depth in both cameras, the same test that picked the true pose from four candidates) and parallax angle (discard the slivers). Calibration $K$ comes from Chapter 12.

def two_view_reconstruct(ptsL, ptsR, K, min_parallax_deg=1.5):
    """Matches + intrinsics -> relative pose and a filtered 3D point set."""
    E, mask = cv2.findEssentialMat(ptsL, ptsR, K,
                                   method=cv2.USAC_MAGSAC,
                                   prob=0.999, threshold=1.0)
    _, R, t, pose_mask = cv2.recoverPose(E, ptsL, ptsR, K, mask=mask)
    keep = pose_mask.ravel().astype(bool)
    pL, pR = ptsL[keep], ptsR[keep]

    P1 = K @ np.hstack([np.eye(3), np.zeros((3, 1))])
    P2 = K @ np.hstack([R, t])
    Xh = cv2.triangulatePoints(P1, P2, pL.T.astype(np.float64),
                                       pR.T.astype(np.float64))
    X = (Xh[:3] / Xh[3]).T

    # cheirality: positive depth in BOTH camera frames
    z1 = X[:, 2]
    z2 = (R @ X.T + t).T[:, 2]
    ok = (z1 > 0) & (z2 > 0)

    # parallax: angle between the two viewing rays at each point
    r1 = X / np.linalg.norm(X, axis=1, keepdims=True)
    C2 = (-R.T @ t).ravel()                  # right camera center
    r2 = X - C2
    r2 /= np.linalg.norm(r2, axis=1, keepdims=True)
    cosang = np.clip((r1 * r2).sum(axis=1), -1, 1)
    ok &= np.degrees(np.arccos(cosang)) > min_parallax_deg

    return R, t, X[ok], ok.sum() / len(X)

R, t, cloud, keep_frac = two_view_reconstruct(ptsL[inl], ptsR[inl], K)
print(f"{len(cloud)} points kept ({100 * keep_frac:.0f}%), "
      f"baseline direction {t.ravel().round(2)}")
# 1147 points kept (83%), baseline direction [-0.97  0.03  0.24]

Two properties of the output deserve flagging. The cloud is scaleless: recoverPose returns $\mathbf{t}$ as a unit vector, so the reconstruction floats at an arbitrary global scale until something metric pins it (a known baseline as in Section 13.5, a known object size, GPS, or an IMU). And it is sparse and view-dependent: points exist only where Chapter 10 found matchable texture. Both limitations are the opening brief for Chapter 14, which chains many pairs, shares scale through overlapping views, and polishes everything with bundle adjustment; the densification story then continues into the multi-view stereo and neural fields of Chapter 27.