Stereo photogrammetry

The hydraulic engineers wanted the depth of the existing drainage ditch where the new canal would cross it, the height of the coulee rim above its floor at three possible crossing sites, and spot elevations on the old canal banks. All of it lay inside the overlap of consecutive corridor photos. Leah set up the first stereo pair on the workstation, put on the viewing glasses, and saw the coulee rise out of the screen in three dimensions. This lesson explains why that works and how it becomes measurement: stereoscopic vision, relative and absolute orientation of a stereo model, the parallax equations of the normal case, the precision of heights and its dependence on the base-to-height ratio, epipolar geometry, and the least-squares space intersection that modern software uses for every point.

Seeing in three dimensions

Human depth perception relies heavily on the slightly different views from the two eyes, about 65 mm apart. The brain fuses the two images, and the difference in the positions of an object in them, the angle of convergence, is perceived as distance. Two aerial photos taken from points hundreds of metres apart are like a pair of eyes separated by the air base, and when each eye sees only its own photo, the brain fuses them into a three-dimensional model of the terrain with greatly exaggerated relief. Viewing systems for photogrammetry deliver the separate images in different ways: mirror stereoscopes for paper prints, anaglyph glasses with red and cyan filters, polarized screens with matching glasses, and active shutter glasses synchronized with the monitor.

The exaggeration is useful. The ratio of the air base to the flying height is typically 0.3 to 0.6 for aerial photography, while the ratio of the eye base to the viewing distance is about 0.15 to 0.25, so the model appears with heights stretched several times relative to horizontal distances. Slopes look steeper than they are, small embankments stand out, and an operator can detect height differences far smaller than could be seen in the field. The same effect means that stereo viewing is a poor guide to actual slope: measurements, not impressions, describe the terrain.

Orienting a stereo model

Before measurements can be made, the two photos must be oriented. In the analogue tradition this was done in two stages, and the stages remain useful for understanding. Relative orientation reproduces the relative position and attitude of the two cameras so that every pair of corresponding rays intersects. It has five independent parameters: the base direction (two angles) and the three relative rotations, the length of the base being arbitrary because it only scales the model. When relative orientation is correct, the y-parallax, the difference in y coordinates of a point in the two photos after rotation into a common frame, is zero everywhere; residual y-parallax shows as an inability to fuse the images.

Absolute orientation then scales, rotates and translates the resulting model into the ground coordinate system with a seven-parameter similarity transformation, three translations, three rotations and a scale, determined from at least two full control points and one height point, in practice more. The two stages together determine twelve parameters, which equal the six exterior orientation parameters of each of the two photos. Modern software usually skips the stages and computes the exterior orientation of all photos directly in aerial triangulation; the stereo workstation simply loads the result.

The normal case

The simplest geometry has both photos exactly vertical, taken from the same altitude, with the x axes parallel to the flight line. This normal case is a useful approximation to real aerial photography and gives closed-form equations. A ground point at elevation \(h\) is imaged at \(x\) on the left photo and \(x'\) on the right, and its x-parallax is

By similar triangles, using the air base \(B\), the flying altitude \(H\) above the datum and the focal length \(f\),

with \(X\) and \(Y\) measured from the left exposure station in the flight direction and across it. These are the parallax equations. Parallax increases with elevation: points closer to the camera are displaced more between the two photos, just as nearby objects shift more than distant ones when you look at them with one eye and then the other. This is the stereo version of relief displacement from Lesson 1.

Heights are often wanted relative to a nearby reference rather than absolutely. Differentiating the parallax equation, a parallax difference \(\Delta p\) between two points corresponds to a height difference

where \(h_r\) is the elevation of the reference point and \(H_g\) the flying height above the ground. A floating mark, two small dots presented one in each image, appears to float above or sink into the terrain as the parallax between the dots is changed; setting it on the ground surface is the operator's way of measuring parallax.

Precision of heights

The parallax equation shows how errors in measured parallax become errors in height. Differentiating \(h = H - Bf/p\) with respect to \(p\) and writing the image scale number \(H_g/f\),

The first factor, the inverse of the base-to-height ratio, expresses the geometry of the intersection: a short base gives rays that meet at a small angle and a poorly defined intersection. The second factor converts image measurement to ground scale. A parallax precision of half a pixel at 10 cm GSD therefore gives a height precision of \(0.05 / (B/H)\) m, about 0.19 m for the corridor flight with 60 % overlap. Horizontal precision, by contrast, depends mainly on the GSD and is better than height precision by roughly the factor \(B/H\). This is why heights from photogrammetry are usually two to four times less precise than horizontal positions, and why lidar is often preferred for bare-earth heights while photogrammetry provides the horizontal detail.

For a given GSD there are two ways to improve heights: a larger base, by reducing the overlap, or a wider-angle camera flown lower, which keeps the GSD and footprint but doubles the base-to-height ratio. Both have costs. Less overlap weakens the block and gives less redundancy for automated matching; a wider angle increases relief displacement and the hidden areas behind buildings and steep slopes. Modern image matching with 80 % overlap compensates for the short base by matching each point in many images rather than two, so that the effective base is the span of all the images that see the point.

A worked program

The program below applies the parallax equations to four points along the corridor, computes the depth of the drainage ditch from a parallax difference, compares height precision for three flight designs, and finally intersects the two rays of one point by least squares.

# Normal-case stereo: heights from parallax, their precision, and a least-squares
# space intersection of one point from two photos
import math

ZREF = 900.0
f = 0.100; B = 520.0                        # focal length (m), air base (m)
ZL = ZREF + 2000.0                          # flying altitude of both exposures
pts = {"canal invert": 0.025976, "canal bank": 0.026000, "coulee rim": 0.026804, "coulee bottom": 0.025490}
print("Point           parallax (mm)   elevation (m)   relative to average terrain (m)")
for name, p in pts.items():
    h = ZL - B * f / p
    print(f"{name:15s} {1000 * p:12.3f}   {h:12.2f}   {h - ZREF:+10.2f}")

# Height difference from a parallax difference (bank vs invert)
pa, pb = pts["canal bank"], pts["canal invert"]
dh = (ZL - B * f / pa) - (ZL - B * f / pb)
print(f"Canal depth from a parallax difference of {1000 * (pa - pb):.3f} mm: {dh:.2f} m")

# Height precision: sigma_Z = (H/B) * (H/f) * sigma_p
print()
print("Height precision for a parallax precision of 0.5 pixel (2.5 microns)")
sp = 2.5e-6
for name, base, ff, H in (("60 % overlap, f = 100 mm", 520.0, 0.100, 2000.0),
                          ("80 % overlap, f = 100 mm", 260.0, 0.100, 2000.0),
                          ("60 % overlap, f = 50 mm", 520.0, 0.050, 1000.0)):
    # same 10 cm GSD and footprint in each case: halving f halves H, the base stays 520 m
    print(f"  {name:26s} H = {H:5.0f} m  B/H = {base / H:.2f}  sigma_Z = {(H / base) * (H / ff) * sp:.3f} m")

# Space intersection: rays from two exposure stations with known orientation (level photos here)
L1 = (10000.0, 20000.0, ZL); L2 = (10000.0 + B, 20000.0, ZL)    # flight line along X
x1, y1 = 0.008257, 0.003327                  # image coordinates on the left photo (m)
x2, y2 = -0.017790, 0.003325                 # on the right photo
# Each ray: P = L + s * (x, y, -f). Least squares for (X, Y, Z, s1, s2) reduces to the
# closest points of two lines; solve the 2 x 2 system for s1, s2.
d1 = (x1, y1, -f); d2 = (x2, y2, -f)
w0 = [a - b for a, b in zip(L1, L2)]
a = sum(c * c for c in d1); b = sum(c * e for c, e in zip(d1, d2)); c = sum(e * e for e in d2)
d = sum(p * q for p, q in zip(d1, w0)); e = sum(p * q for p, q in zip(d2, w0))
s1 = (b * e - c * d) / (a * c - b * b); s2 = (a * e - b * d) / (a * c - b * b)
P1 = [l + s1 * k for l, k in zip(L1, d1)]; P2 = [l + s2 * k for l, k in zip(L2, d2)]
P = [(u + v) / 2 for u, v in zip(P1, P2)]
gap = math.dist(P1, P2)
print()
print(f"Intersected point X = {P[0]:.3f}  Y = {P[1]:.3f}  Z = {P[2]:.3f}  (rays miss by {100 * gap:.1f} cm)")
print(f"Parallax {1000 * (x1 - x2):.3f} mm gives Z = {ZL - B * f / (x1 - x2):.3f}; y-parallax {1e6 * (y1 - y2):+.1f} um")

Point           parallax (mm)   elevation (m)   relative to average terrain (m)
canal invert          25.976         898.15        -1.85
canal bank            26.000         900.00        +0.00
coulee rim            26.804         959.99       +59.99
coulee bottom         25.490         859.98       -40.02
Canal depth from a parallax difference of 0.024 mm: 1.85 m

Height precision for a parallax precision of 0.5 pixel (2.5 microns)
  60 % overlap, f = 100 mm   H =  2000 m  B/H = 0.26  sigma_Z = 0.192 m
  80 % overlap, f = 100 mm   H =  2000 m  B/H = 0.13  sigma_Z = 0.385 m
  60 % overlap, f = 50 mm    H =  1000 m  B/H = 0.52  sigma_Z = 0.096 m

Intersected point X = 10164.842  Y = 20066.400  Z = 903.609  (rays miss by 4.0 cm)
Parallax 26.047 mm gives Z = 903.609; y-parallax +2.0 um

The parallax heights reproduce the terrain: the canal invert at 898.15 m, the bank at 900.00 m, the rim at 959.99 m and the coulee bottom at 859.98 m. The ditch depth of 1.85 m comes from a parallax difference of only 24 µm, about five pixels; an error of half a pixel in that difference would change the depth by about 19 cm, the precision computed for this geometry. The intersected point has X = 10 164.842, Y = 20 066.400 and Z = 903.609 m. Its two rays miss each other by 4 cm, an amount consistent with image errors of a micrometre or two, and the residual y-parallax is +2.0 µm.

Space intersection

In general, photos are not in the normal case, and the ground point is found by space intersection: the exterior orientation of two or more photos is known, the point is measured in each, and its three ground coordinates are the unknowns. Each photo gives two collinearity equations, so two photos give four equations in three unknowns, one degree of freedom, which appears as the y-parallax. With more photos the redundancy grows by two per photo. The solution is by iterated least squares, or for two rays by finding the closest points of the two lines and taking their midpoint, as in the program. The residual of the intersection is a check on the measurement and on the orientation: large residuals mean a mismeasured point or a wrongly oriented photo.

Epipolar geometry

The epipolar plane of a point contains the two perspective centres and the point; it cuts each image along an epipolar line. A point's image in the right photo must lie on the epipolar line defined by its image in the left photo. This reduces the search for corresponding points from two dimensions to one and is used by every matching algorithm. Stereo workstations resample both images so that epipolar lines become horizontal rows, producing normalized or epipolar images in which the normal-case equations apply exactly and the operator, or the matching algorithm, need only move along rows. All epipolar lines in an image pass through the epipole, the image of the other camera's perspective centre; for near-vertical aerial photos the epipoles lie far outside the format and the lines are nearly parallel.

Stereo compilation

Stereo compilation is the measurement of features in a stereo model: break lines along the top and toe of the canal banks, the edges of roads and buildings, spot heights on the coulee floor, and contours traced by keeping the floating mark on the ground at a constant height. Although automated matching now generates most surface models, human stereo compilation is still used for break lines, for features that matching handles poorly such as narrow ditches and bridge decks, and for quality control. Leah compiled the ditch crossing by hand: the matched surface smoothed the 1.85 m ditch into a gentle swale, while her break lines along its top and bottom preserved its shape for the hydraulic model.