Photogrammetry

Photogrammetry is measurement from photographs: the coordinates of points, lines and surfaces are computed from their images, in a stated reference system and with a stated accuracy. It is older than computer vision and grew up with aerial mapping, and its central computation, the bundle adjustment, is the one that structure from motion adopted and later recognized (see computer vision). Photogrammetry is called aerial when the images are taken from aircraft looking down at the ground, and close-range (or terrestrial) when the camera is near the object, as for a façade, a rock face, an excavation or an industrial part. Workflows that run structure from motion and multi-view stereo on ordinary photographs are called SfM photogrammetry or SfM-MVS in the surveying and geoscience literature. This article treats photogrammetry as a surveying application: how image networks are tied to the ground by control points, how the camera is calibrated within the adjustment, how accuracy is defined and tested at independent check points, which products result, and where the method is used. The reconstruction algorithms themselves belong to structure from motion, bundle adjustment and multi-view stereo, and the planning and capture of drone imagery, with ground sampling distance and image overlap, to aerial and drone image analysis.

Problem

A photogrammetric survey answers a measurement question: where is this point, how high is this surface, how much material has moved, and how sure are we? Images alone determine a scene only up to a similarity transformation, so a reconstruction from them has a shape but no position, orientation or scale on the Earth (Szeliski, 2022, §11.3). Knowledge of absolute positions or distances, which photogrammetry calls ground control points, is required to fix them, and Szeliski adds that where absolute metric accuracy is required, as in photogrammetry, the cameras should be pre-calibrated and ground control points used to pin down the reconstruction, whereas visualization can rely on a simple camera model (Szeliski, 2022, §11.3).

Two kinds of quality follow from this, and the literature keeps them apart. Internal precision describes how well the image network fixes the relative shape; it can be read from the residuals and covariances of the adjustment. Accuracy describes how close the computed coordinates are to the truth in the target reference system, and can only be measured against independent measurements of higher accuracy. Sanz-Ablanedo, Chandler, Rodríguez-Pérez and Ordóñez (2018) point out that residuals of the bundle adjustment after the model is rotated and scaled are internal precision rather than accuracy, and they do not use them as an accuracy measure. James, Chandler, Eltner and co-authors (2019) separate a second pair: bias, the systematic part of the error (for example the mean error), and precision, the random part (for example its standard deviation).

The geometry of the image network matters as much as the camera. Remondino, Barazzetti, Nex, Scaioni and Sarazzi (2011) contrast the two settings. In aerial photogrammetry, the regular blocks of vertical images are oriented by automatic aerial triangulation. In close-range work, automation is harder because the image scale is large and often varying, the geometry is convergent, the overlap irregular and the changes in geometry and radiometry strong; they add that blocks from unmanned aircraft are in many cases more similar to close-range than to aerial blocks.

History in outline

Photogrammetry supplied several of the methods that computer vision later took up. Szeliski (2022) records that camera calibration was first studied in photogrammetry (§11.6); that the earliest stereo matching algorithms were developed there to build topographic elevation maps from overlapping aerial images, after a period in which operators traced constant-elevation contours on stereo plotters (introduction to Chapter 12); and that “bundle adjustment” is the photogrammetric name for the nonlinear minimization of reprojection error (§11.4.2). The bundle adjustment and lens distortion articles trace the bundle method and the Brown–Conrady distortion model to photogrammetry. Szeliski also notes that in classic photogrammetry the reporting of precision is almost always considered mandatory (Appendix B.6), while many computer vision applications, by his account, ignore uncertainty modelling.

The recent change came from the other direction. James et al. (2019) describe structure-from-motion photogrammetry as combining digital photogrammetry with the flexibility and ease of use of multi-view computer vision, list among its novelties automatic orientation of image networks from aerial or terrestrial platforms and automated camera calibration suited to consumer cameras, in low-cost and sometimes free software, and observe that most of its new users in geomorphology have little or no formal training in photogrammetry. Their guidelines for the journal Earth Surface Processes and Landforms exist to bring the field’s error analysis to those users.

Visual Data

Images. Aerial blocks are strips of near-vertical (nadir) images with overlap along and between the strips, from manned aircraft or, at low altitude, from drones; how such blocks are planned is described in aerial and drone image analysis. Oblique images, taken with the camera tilted, add views of façades and steep slopes and strengthen the geometry. Close-range networks surround or face the object with convergent images. James et al. (2019) count metric cameras among the expensive technology that kept photogrammetry a specialist method, and contrast them with the off-the-shelf consumer cameras whose calibration SfM photogrammetry automates. A consumer camera can change its geometry between uses, which is why Sanz-Ablanedo et al. (2018) gave each of their six flights its own camera model, judging that a landing impact might have changed the camera.

Ground control. A ground control point (GCP) is a point recognizable in the images whose coordinates are measured on the ground, typically with survey-grade satellite positioning or a total station. Giordan et al. (2020) distinguish natural points (corners, road markings and other features that are well defined as points rather than blobs), artificial markers placed on the ground, whose centre is precisely defined, and coded markers whose pattern also identifies them. The same kind of point serves as a check point when its coordinates are withheld from the adjustment and used only to test it. James et al. (2019) list scale bars, as well as ground control points and camera positions or orientations, among the control measurements used to scale and georeference results.

Camera positions. Satellite positioning on the aircraft gives a position for each image. Sanz-Ablanedo et al. (2018) explain why this is hard to use as control: the camera moves at several metres per second, so triggering must be synchronized with the receiver, and no repeated measurements at the same point are possible. They conclude from earlier studies that direct georeferencing achieves only decimetre to metre accuracies, and report that receivers with real-time kinematic corrections came closer to control-point results horizontally but less so vertically; one of the studies they cite concluded that ground control should be used when a vertical RMSE better than 12 cm is needed.

Vision Tasks Involved

Photogrammetry and computer vision name the same computations differently, and reading either literature needs both sets of names. The table is our mapping; the five-point algorithm, perspective-n-point and Kabsch algorithm articles give the photogrammetric names relative orientation, space resection and absolute orientation themselves:

Photogrammetry	Computer vision
Interior orientation, camera calibration	Intrinsics and lens distortion (camera calibration)
Exterior orientation	Camera pose
Relative orientation	Relative pose from the essential matrix (five-point algorithm)
Space resection	Absolute pose (perspective-n-point)
Space intersection	Triangulation
Absolute orientation	Similarity alignment of point sets (Kabsch algorithm)
(Aerial) triangulation, bundle block adjustment	Bundle adjustment
Free network, inner-constraints adjustment	Bundle adjustment with a gauge fixed by constraints
Tie points	Feature tracks
Dense image matching	Multi-view stereo

The vision tasks are those of structure from motion and multi-view stereo. What photogrammetry adds is the step that ties the network to the ground and the statistics that say how well it does so.

Common Methods

Georeferencing

Sanz-Ablanedo et al. (2018) define georeferencing as transforming the arbitrary datum of an initial SfM model into a predefined coordinate reference system. Direct georeferencing uses the measured positions (and possibly orientations) of the cameras; indirect georeferencing uses ground control points identified in the images. Ground control can enter in two ways.

Similarity after adjustment. The model is reconstructed from the images alone and then moved, rotated and scaled onto the control coordinates by a seven-parameter similarity transformation, which Sanz-Ablanedo et al. (2018) call a conformal transformation:

$$\mathbf{X}_{\text{world}} = s\,R\,\mathbf{X}_{\text{model}} + \mathbf{t}.$$

Given three or more non-collinear points known in both frames, the least-squares $s$, $R$ and $\mathbf{t}$ have the closed form of the Kabsch algorithm with Umeyama’s scale, the problem photogrammetry calls absolute orientation. COLMAP’s geo-registration estimates such a similarity: given the 3D positions of the camera centres of at least three registered images, it estimates the similarity between model and target frame with RANSAC, using a user-supplied error threshold to reject outliers (COLMAP 4.2.1 documentation). A similarity keeps the model’s shape. Sanz-Ablanedo et al. (2018) therefore note that it cannot fully correct deformations of the model, since it only translates, rotates and scales it to minimize the residuals at the control points.

Control inside the adjustment. The second approach lets the control reshape the model: the surveyed coordinates become observations in the bundle adjustment, weighted by their precision, alongside the image observations. With image observations $\mathbf{x}_{ij}$ of precision $\sigma_{\text{img}}$ and control points $g$ with surveyed coordinates $\hat{\mathbf{X}}_g$ of precision $\sigma_g$ per axis, the adjustment minimizes

$$E = \sum_{(i,j)} \frac{\big\lVert \mathbf{x}_{ij} - \pi(\mathbf{K}, R_i, \mathbf{C}_i, \mathbf{X}_j) \big\rVert^2}{\sigma_{\text{img}}^2} \;+\; \sum_{g} \frac{\big\lVert \mathbf{X}_g - \hat{\mathbf{X}}_g \big\rVert^2}{\sigma_g^2},$$

over the cameras $(R_i, \mathbf{C}_i)$, the points $\mathbf{X}_j$ (tie, control and check points alike) and, if the camera is self-calibrated, its intrinsics $\mathbf{K}$. Under independent Gaussian errors with these standard deviations this is the maximum-likelihood estimate, as in the bundle adjustment article. The control terms also remove the seven-parameter gauge freedom that images alone leave (our reading of the gauge discussion in bundle adjustment): three or more well-placed control points fix position, orientation and scale. COLMAP’s pose_prior_mapper applies the same idea to camera positions from the images’ GPS metadata, as prior position constraints during incremental mapping; the user can override the priors’ covariance with standard deviations per axis that default to 1 m (COLMAP 4.2.1 documentation).

James, Robson, d’Oleire-Oltmanns and Niethammer (2017) explain why the weights matter. Control measurements define the datum, and, if well distributed, constrain the systematic errors in shape that accumulate across extended networks. If the control precision is overstated, the shape is distorted by overfitting to errors in the control at the expense of the tie points; if it is understated, control does less to suppress systematic error. They add that the weights of the image observations of tie points and of control points can matter as much as the precision assigned to the field survey. In the software they used (PhotoScan 1.1.6), the defaults described image observations of control points as precise to 0.1 pixel and tie points to 4.0 pixels, with 5 mm for the control coordinates. On their erosion-gully survey in Morocco these defaults gave a ratio of about 6 between the errors on check points and on control points, which they read as overweighting of the control; values taken from the actual image residuals of exploratory adjustments, 0.69 pixel for tie points and 0.50 pixel for control observations, brought the ratio to about 1.3, and the refined settings removed step-like artefacts of up to about 50 mm that had followed changes in image overlap.

Free network first. The same study screened the data in a sequence they draw from conventional photogrammetric quality checks: an adjustment with tie points only, then one including the image observations of the GCPs without their ground coordinates (a free or inner-constraints adjustment, in which the GCP image measurements can be checked against the network), and only then one with the ground coordinates as control, whose discrepancies expose misidentified or moved control points (James et al., 2017). One of their 30 GCPs, with acceptable image residuals, introduced a systematic tilt of about 20 mm across the model and was removed; on a landslide survey with 291 targets, 47 were removed with a threshold just above four times the assumed GCP precision of 0.05 m, probably because parts of the slide had moved differently during the several days of image capture.

Camera calibration and self-calibration

The camera model is the pinhole camera with focal length (photogrammetry’s principal distance), principal point and lens distortion, sometimes with affinity and skew terms. It can be fixed from a separate calibration or estimated within the adjustment, which is self-calibration. Remondino et al. (2011) write that photogrammetric applications prefer to separate calibration from orientation, because the two need different network geometries, also for drone blocks without cross strips, and that self-calibrating bundle adjustment is used where lower accuracy suffices. SfM photogrammetry usually self-calibrates, often with many parameters: Sanz-Ablanedo et al. (2018) estimated focal length, principal point, four radial and four tangential distortion coefficients and two affinity and skew coefficients in every adjustment, and report correlations as high as −0.98 and −0.99 between successive radial coefficients.

Self-calibration can overfit. James et al. (2017) compared four sets of camera parameters inside a Monte Carlo loop that reran the adjustment with different random halves of the GCPs as control: (A) focal length only; (B) adding the principal point and three radial terms; (C) adding two tangential terms; (D) four radial and two tangential terms, aspect ratio and skew. A parameter set was accepted only if it lowered the RMSE on both control and check points without substantially raising the ratio between them. On the gully survey, errors fell from A to C, but D raised the ratio of check to control error, the signature of overfitting the control. They also cite laboratory calibrations in which only the principal point offsets and two radial parameters of DSLRs with 28 mm lenses exceeded their own precision estimates. Their conclusion on the trade-off: with dense ground control, self-calibrating surveys are likely to outperform pre-calibrated ones, while surveys with few GCPs should use pre-calibrated cameras, whose calibration must then be better than the precision the network needs, which is hard for consumer cameras whose internal geometry is not stable.

Doming. The best documented systematic error of SfM photogrammetry is a dome- or bowl-shaped deformation of the surface model. Sanz-Ablanedo et al. (2018) attribute it to projective coupling or correlation between the inner and exterior orientation parameters, which leads to an inaccurate lens model and a deformed model, and note that processing vertical imagery with automated SfM procedures is recognized to produce an error surface in the form of a systematic dome; they cite several studies, among them James and Robson (2014), which could not be read for this article. In their own block, five control points placed semi-randomly towards the edge of a 1,225 ha site gave an RMSE of 0.03 m at the control points while many check points were more than 1 m off in height, with the dome’s maximum as far as possible from all the control points and not coinciding with the mine pit in the middle of the site. James et al. (2017) saw doming on their gully survey when the control was given very little weight, and saw it return with their most complex camera model. Two remedies recur in these sources. The first is more and better spread control; the second, which James et al. (2017) attribute to James and Robson (2014), is to add oblique images, which strengthen the network geometry and reduce errors in the estimated camera parameters, minimizing systematic deformation even without control. OpenDroneMap documents a third: its option use-fixed-camera-params turns off camera parameter optimization during bundle adjustment, which its documentation describes as sometimes useful for results that exhibit doming or bowling, or for rolling-shutter cameras (OpenDroneMap 3.6.0 documentation).

Products

After orientation, dense matching produces the surface products. Remondino et al. (2011) describe automated matching delivering a dense point cloud of the surface, the digital surface model, which is interpolated, simplified and textured, followed by orthoimages; for man-made features and vector data, interactive measurement was still the norm in their account. The products are:

• Point clouds of the surface, georeferenced, often with colour (point cloud).
• Digital surface models (DSMs), the height of the visible surface including buildings and vegetation, and digital terrain models (DTMs) of the bare ground, usually as height grids (a depth map in map geometry).
• Orthophotos, images resampled onto the surface model so that every pixel has map geometry; aerial and drone image analysis explains why they differ from homography mosaics and how they fail at building edges.
• Textured meshes for visualization and documentation (polygon mesh).

The dense steps are those of multi-view stereo and stereo matching, with semi-global matching among the methods Remondino et al. (2011) list as state of the art for multi-image matching. James et al. (2019) warn that analyses of these products must keep the uncertainty of the photogrammetric processing, including possible systematic error, and consider the smoothing and filtering that dense matching applies when they derive resolutions, roughness or surface change.

Evaluation

Check points and the NSSDA statistics

The United States National Standard for Spatial Data Accuracy (NSSDA, Federal Geographic Data Committee, 1998) is a readable official example of how positional accuracy is tested. It applies to georeferenced maps and digital data derived from aerial photographs, satellite imagery and ground surveys, uses the root-mean-square error to estimate positional accuracy, and reports accuracy in ground distances at the 95% confidence level. For $n$ check points with dataset coordinates $(x_i, y_i, z_i)$ and coordinates $(\hat x_i, \hat y_i, \hat z_i)$ from an independent source of higher accuracy,

$$\mathrm{RMSE}_x = \sqrt{\tfrac{1}{n}\textstyle\sum_i (x_i - \hat x_i)^2}, \qquad \mathrm{RMSE}_r = \sqrt{\mathrm{RMSE}_x^2 + \mathrm{RMSE}_y^2}, \qquad \mathrm{RMSE}_z = \sqrt{\tfrac{1}{n}\textstyle\sum_i (z_i - \hat z_i)^2},$$

with $\mathrm{RMSE}_y$ defined like $\mathrm{RMSE}_x$. Assuming that systematic errors have been removed as far as possible and that the errors are normally distributed and independent in $x$ and $y$, the standard’s horizontal accuracy at 95% confidence is $1.7308\,\mathrm{RMSE}_r$ when $\mathrm{RMSE}_x = \mathrm{RMSE}_y$ (from a factor 2.4477 for the circular error), with an approximation $2.4477 \times 0.5(\mathrm{RMSE}_x + \mathrm{RMSE}_y)$ when the ratio of the smaller to the larger component lies between 0.6 and 1.0; vertical accuracy is $1.9600\,\mathrm{RMSE}_z$ (Appendix 3-A). The check points must be well-defined points, acquired separately from the data used in the aerotriangulation; at least 20 must be tested, distributed to reflect the area and the distribution of error, so that at 95% confidence one of 20 may exceed a product’s threshold. For a rectangular area believed to have uniform accuracy, the standard’s informative guideline spaces check points at least 10% of the diagonal apart with at least 20% in each quadrant (Appendix 3-C). The NSSDA sets no thresholds: agencies and users decide what accuracy their application needs. The NSSDA was written to update the US National Map Accuracy Standards of 1947, and it took the ASPRS Accuracy Standards for Large-Scale Maps of 1990 as its basis.

What RMSE hides

RMSE at check points is the common currency, and both sources read here qualify it. James et al. (2019) ask authors to avoid statistics that mix systematic and random error, such as RMSE alone, and to report a measure of bias (such as the mean error) and of precision (such as its standard deviation), to examine how the error varies over the area, and to non-dimensionalize results by the mean observation distance, the survey extent or the ground sampling distance so that surveys can be compared. Where no check data exist, they propose a split test: two datasets of the same area from different survey designs or dates, compared over zones known to be stable. James et al. (2017) showed the limit of RMSE directly: on the gully survey, adjustments with four to fifteen well-distributed control points gave check-point RMSEs of 47–50 mm with fixed and with self-calibrated cameras alike, yet comparing their surface models with one adjusted to all GCPs revealed systematic differences, linear slopes and doming, that the RMSE values concealed, most apparent with self-calibration.

Control points are not check points. Sanz-Ablanedo et al. (2018) found that the RMSE at control points stays far below that at check points, because the model adapts to whatever control it is given. With three to five control points the control RMSE was sometimes a few millimetres while the check-point RMSE exceeded 1 m and reached 8 m. If no check points are available, they suggest estimating the real 3D accuracy by multiplying the control RMSE by an overestimation factor, about 3 for projects with many GCPs (more than 3.5 per 100 photos in their study) and 4 to 8 for projects with few (fewer than 2 per 100), and they note that the factor is larger for vertical and smaller for horizontal errors. James et al. (2019) state the same rule as a reporting requirement: data quality must not be judged from the error at control measurements alone.

Precision is not accuracy. On an archaeological block at Veio, Italy, Remondino et al. (2011) report that the covariance matrix of the bundle solution gave a theoretical precision of the object coordinates of ±0.6 cm horizontally and ±2.3 cm vertically, while the RMSE on five targets was 4, 3 and 7 cm in X, Y and Z, which they attribute to the manual measurement of target centres that were not clearly visible in the images.

How many control points, and where

Sanz-Ablanedo et al. (2018) ran the largest systematic test among the sources read here. On a 1,225 ha open-cast coal mining area in León, Spain, with 558 m of relief, 2,514 images from a fixed-wing drone at an average ground sampling distance (GSD) of 6.86 cm, and 102 targets measured by RTK satellite positioning (the caption of their Figure 7 speaks of 2,415 images), they ran 3,465 self-calibrating bundle adjustments: 35 random selections of control points for each number from 3 to 101, the rest serving as check points. Their main findings:

• With few control points (10–20), the 3D check-point RMSE was over 31 cm, about five times the average GSD; with 50–60 it was about 16 cm (3 GSD); with 90–100 it converged slowly to about 12 cm, about twice the GSD, and did not improve further.
• Vertical errors at check points were about 2.5 times the easting or northing errors. Horizontal accuracy reached the GSD with about 2.5–3 control points per 100 photos and improved little beyond; vertical accuracy kept improving towards 1.5 times the GSD, reached with 4 per 100 photos. Below about one control point per 100 photos accuracy degraded rapidly.
• For a given number of points, well-distributed control gave about half the RMSE of poorly distributed control. Concentrating points in parts of the area, leaving gaps, or placing them only on the periphery or only in the centre all performed poorly; the authors recommend an even distribution, ideally a triangular grid, which minimizes the largest distance from any point to its nearest control point.
• Within the block, points with a coarser GSD (lower ground, farther from the camera) were less accurate, but by about half the expected amount, which the authors suggest may be because such points appear in more images.

James et al. (2017) reached a complementary conclusion with Monte Carlo resampling of their two surveys: once an appropriate camera model and weights were used, roughly half the GCPs could have been deployed with little loss, and the GCP coverage should extend to the extremities of the area so that surfaces are interpolated within the control rather than extrapolated beyond it. The two studies differ in setting, a 1,225 ha site with large relief and 102 targets against a gully survey with 30 targets and a landslide with 291, and their numbers are specific to those sites, as Sanz-Ablanedo et al. say of theirs.

Example: Control Points and Check Points in a Simulated Block

The script simulates a survey to show how the number and placement of control points change what check points measure. It is a simulation for understanding accuracy, not advice on how to fly or survey. Sixty nadir images, 3000 × 2000 pixels with a focal length of 2,500 pixels, are taken from 100 m above a datum on six parallel lines over synthetic terrain with hills up to about 12 m, which gives a GSD of 4 cm at the datum. The true camera has radial distortion ($k_1 = -0.08$, $k_2 = 0.02$ on normalized coordinates) and a principal point a few pixels off centre. About 2,000 terrain points are projected into every image that sees them; 30 of them, on a jittered grid, are check points, and five sets of ground control points are compared: four corners, the corners plus the centre, a 3 × 3 and a 4 × 4 spread grid, and four points clustered in one corner. Image observations of all points receive Gaussian noise of 0.5 pixel per coordinate, and surveyed control coordinates 1 cm per axis.

The adjustment minimizes the weighted cost $E$ above with the true standard deviations as weights, over all camera poses, all points and, under self-calibration, a shared focal length, principal point, $k_1$ and $k_2$. It starts from camera centres perturbed by 2 m, angles by 1° and points by 1 m, standing in for an SfM model already brought near the ground frame. The solver is a compact Levenberg–Marquardt with an analytic Jacobian and the Schur complement over the 3 × 3 point blocks described in bundle adjustment. Conventions: the world frame has $z$ up; $R_i$ maps world to camera coordinates, $\mathbf{p} = R_i(\mathbf{X} - \mathbf{C}_i)$, with the camera looking down its $+z$ axis; rotation updates are applied on the left, $R_i \leftarrow \exp([\boldsymbol{\delta}]_\times) R_i$ (the bundle adjustment article writes the update on the right; either parametrizes the same rotations). Accuracy is scored as the NSSDA defines it, by $\mathrm{RMSE}_r$ and $\mathrm{RMSE}_z$ at the check points, together with the signed mean vertical error (bias, as James et al. (2019) ask) and the RMSE at the control points themselves. Two calibration choices are printed, because the sources leave the choice to the survey: a camera calibrated beforehand and held fixed at its true values, and self-calibration. Each configuration runs with noise seeds 0–4, and the script prints the mean and the range over them. Before the first run, a pilot at these sizes confirmed that every configuration recovers the check points exactly from noise-free data under both choices, a result the script prints again for one case; the plan’s first solver, SciPy’s least_squares, was replaced before the first run by the one shown because it needed hundreds of evaluations per adjustment on noisy data. The plan expected that under self-calibration “corners-only and clustered show larger vertical error with a dome-like pattern”; the vertical error grids for those two configurations were added to the printout before the first run, but after a pilot’s seed-0 accuracy figures for them had been seen. The script shown is the first run unchanged. It takes from about 20 seconds to about a minute with NumPy 2.3.5 and SciPy 1.16.3, depending on the machine, and repeated runs gave identical output.

"""Ground control and check points in a simulated nadir image block.
Bundle adjustment with GCPs as weighted observations; accuracy measured at
independent check points (NSSDA: RMSE_r horizontal, RMSE_z vertical)."""
import numpy as np
import scipy.sparse as sp
from scipy.spatial.transform import Rotation

W, H = 3000, 2000
K_TRUE = np.array([2500.0, 1503.0, 998.0, -0.08, 0.02])   # f, cx, cy, k1, k2 (pixels; k on normalized coords)
ALT = 100.0                                                # camera height above z = 0 (m)
SIG_IMG, SIG_GCP = 0.5, 0.01                               # image noise (px), survey noise (m)
R0 = np.diag([1.0, -1.0, -1.0])                            # nadir view: world -> camera

CONFIGS = {
    "4 corners": [(15, 15), (285, 15), (15, 185), (285, 185)],
    "5 corners+centre": [(15, 15), (285, 15), (15, 185), (285, 185), (150, 100)],
    "9 spread (3x3)": [(x, y) for x in (15, 150, 285) for y in (15, 100, 185)],
    "16 spread (4x4)": [(x, y) for x in (15, 105, 195, 285) for y in (15, 72, 128, 185)],
    "4 clustered": [(15, 15), (75, 15), (15, 60), (75, 60)],
}


def terrain(x, y):
    hills = [(80, 60, 40, 8.0), (220, 150, 50, 12.0), (260, 40, 30, 6.0), (120, 170, 35, 5.0)]
    return sum(h * np.exp(-((x - a) ** 2 + (y - b) ** 2) / (2 * s * s)) for a, b, s, h in hills)


def build_scene():
    rng = np.random.default_rng(12345)                     # fixed layout seed
    R, C = [], []
    for x in np.arange(0, 301, 60.0):                      # 6 lines along y
        for y in np.arange(-12, 205, 24.0):                # 10 images per line
            C.append(np.array([x, y, ALT]) + rng.normal(0, 0.5, 3))
            R.append(Rotation.from_rotvec(np.deg2rad(rng.normal(0, 1.0, 3))).as_matrix() @ R0)
    tie = np.c_[rng.uniform(-40, 340, 2000), rng.uniform(-40, 240, 2000)]
    gx, gy = np.meshgrid(37.5 + 45 * np.arange(6), 30 + 35 * np.arange(5))
    chk = np.c_[gx.ravel(), gy.ravel()] + rng.uniform(-5, 5, (30, 2))
    return np.array(R), np.array(C), tie, chk


def project(K, R, C, X, jac=False):
    """Pixels of points X[k] in cameras (R[k], C[k]); optional Jacobians."""
    p = np.einsum("nij,nj->ni", R, X - C)
    a, b, c = p.T
    xn, yn = a / c, b / c
    r2 = xn * xn + yn * yn
    g = 1 + K[3] * r2 + K[4] * r2 * r2
    uv = np.c_[K[0] * g * xn + K[1], K[0] * g * yn + K[2]]
    if not jac:
        return uv, c
    dg = 2 * (K[3] + 2 * K[4] * r2)                        # dg/dxn = dg*xn, dg/dyn = dg*yn
    D = np.empty((len(a), 2, 2))                           # d(u,v)/d(xn,yn)
    D[:, 0, 0] = K[0] * (g + dg * xn * xn)
    D[:, 0, 1] = K[0] * dg * xn * yn
    D[:, 1, 0] = D[:, 0, 1]
    D[:, 1, 1] = K[0] * (g + dg * yn * yn)
    P = np.zeros((len(a), 2, 3))                           # d(xn,yn)/dp
    P[:, 0, 0] = 1 / c
    P[:, 1, 1] = 1 / c
    P[:, 0, 2] = -a / c ** 2
    P[:, 1, 2] = -b / c ** 2
    DP = D @ P
    J_X = DP @ R                                           # d uv / d X
    px = np.zeros((len(a), 3, 3))                          # -[p]x : dp/d(rotation increment)
    px[:, 0, 1], px[:, 0, 2], px[:, 1, 2] = c, -b, a
    px[:, 1, 0], px[:, 2, 0], px[:, 2, 1] = -c, b, -a
    J_cam = np.concatenate([DP @ px, -J_X], axis=2)        # (rotation, centre)
    J_K = np.stack([np.c_[g * xn, np.ones_like(a), np.zeros_like(a), K[0] * xn * r2, K[0] * xn * r2 ** 2],
                    np.c_[g * yn, np.zeros_like(a), np.ones_like(a), K[0] * yn * r2, K[0] * yn * r2 ** 2]], 1)
    return uv, c, J_cam, J_K, J_X


def adjust(R, C, K, X, obs, uv_obs, gcp_ids, gcp_obs, selfcal, iters=100):
    """Levenberg-Marquardt over camera rotations/centres, shared intrinsics (if
    self-calibrating) and points; GCPs enter as weighted coordinate observations."""
    m, n, nk = len(R), len(X), (5 if selfcal else 0)
    ci, pj = obs[:, 0], obs[:, 1]
    na = 6 * m + nk
    rows = np.arange(2 * len(obs)).reshape(-1, 2)

    def resid(R, C, K, X):
        uv, _ = project(K, R[ci], C[ci], X[pj])
        return np.r_[((uv - uv_obs) / SIG_IMG).ravel(), ((X[gcp_ids] - gcp_obs) / SIG_GCP).ravel()]

    def jacobian(R, C, K, X):
        _, _, Jc, Jk, Jx = project(K, R[ci], C[ci], X[pj], jac=True)
        r_, c_, v_ = [], [], []
        for blk, col0 in ((Jc, 6 * ci), (Jx, na + 3 * pj)):
            w = blk.shape[2]
            r_.append(np.repeat(rows[:, :, None], w, 2).ravel())
            c_.append((col0[:, None, None] + np.arange(w)[None, None, :]).repeat(2, 1).ravel())
            v_.append((blk / SIG_IMG).ravel())
        if selfcal:
            r_.append(np.repeat(rows[:, :, None], nk, 2).ravel())
            c_.append(np.broadcast_to(6 * m + np.arange(nk), (len(obs), 2, nk)).ravel())
            v_.append((Jk / SIG_IMG).ravel())
        g_rows = 2 * len(obs) + np.arange(3 * len(gcp_ids))
        r_.append(g_rows)
        c_.append((na + 3 * gcp_ids[:, None] + np.arange(3)).ravel())
        v_.append(np.full(len(g_rows), 1 / SIG_GCP))
        return sp.csr_matrix((np.concatenate(v_), (np.concatenate(r_), np.concatenate(c_))),
                             shape=(2 * len(obs) + 3 * len(gcp_ids), na + 3 * n))

    def update(R, C, K, X, d):
        dc = d[:6 * m].reshape(m, 6)
        R2 = Rotation.from_rotvec(dc[:, :3]).as_matrix() @ R
        K2 = K + d[6 * m:na] if selfcal else K
        return R2, C + dc[:, 3:], K2, X + d[na:].reshape(n, 3)

    r = resid(R, C, K, X)
    cost, lam = 0.5 * r @ r, 1e-3
    for it in range(iters):
        J = jacobian(R, C, K, X)
        Hm = (J.T @ J).tocsr()
        g = -(J.T @ r)
        Haa = Hm[:na, :na].toarray()
        Hab = Hm[:na, na:]
        Vb = Hm[na:, na:].tobsr(blocksize=(3, 3))                        # point blocks (block diagonal)
        Vb = Vb.data[np.argsort(np.repeat(np.arange(n), np.diff(Vb.indptr)) * n + Vb.indices)]
        while True:
            Ua = Haa + lam * np.diag(np.diag(Haa))
            Vd = Vb + lam * Vb * np.eye(3)[None]
            Vi = np.linalg.inv(Vd)
            Vinv = sp.bsr_matrix((Vi, np.arange(n), np.arange(n + 1)), shape=(3 * n, 3 * n)).tocsr()
            S = Ua - (Hab @ Vinv @ Hab.T).toarray()                    # reduced camera system
            da = np.linalg.solve(S, g[:na] - Hab @ (Vinv @ g[na:]))
            db = Vinv @ (g[na:] - Hab.T @ da)                           # back-substitution
            cand = update(R, C, K, X, np.r_[da, db])
            r2 = resid(*cand)
            c2 = 0.5 * r2 @ r2
            if c2 < cost:
                break
            lam *= 10
            if lam > 1e12:
                return R, C, K, X, cost, it
        R, C, K, X = cand
        rel = (cost - c2) / cost
        r, cost, lam = r2, c2, max(lam / 10, 1e-12)
        if rel < 1e-12:
            break
    return R, C, K, X, cost, it + 1


def experiment(config, selfcal, seed, noise=True):
    R_t, C_t, tie_xy, chk_xy = build_scene()
    gcp_xy = np.array(CONFIGS[config], float)
    xy = np.r_[tie_xy, chk_xy, gcp_xy]
    X_t = np.c_[xy, terrain(xy[:, 0], xy[:, 1])]
    n_tie, n_chk, n_gcp = len(tie_xy), len(chk_xy), len(gcp_xy)
    obs = []
    for i in range(len(R_t)):                                          # visibility in the true images
        uv, z = project(K_TRUE, np.repeat(R_t[i][None], len(X_t), 0), np.repeat(C_t[i][None], len(X_t), 0), X_t)
        ok = (z > 0) & (uv[:, 0] >= 0) & (uv[:, 0] < W) & (uv[:, 1] >= 0) & (uv[:, 1] < H)
        obs += [(i, j) for j in np.flatnonzero(ok)]
    obs = np.array(obs)
    cnt = np.bincount(obs[:, 1], minlength=len(X_t))
    assert (cnt[n_tie:] >= 2).all()
    keep = cnt >= 2
    remap = np.cumsum(keep) - 1
    X_t, obs = X_t[keep], obs[keep[obs[:, 1]]]
    obs[:, 1] = remap[obs[:, 1]]
    n = len(X_t)
    chk_ids, gcp_ids = np.arange(n - n_chk - n_gcp, n - n_gcp), np.arange(n - n_gcp, n)

    rng = np.random.default_rng(seed)
    uv_t, _ = project(K_TRUE, R_t[obs[:, 0]], C_t[obs[:, 0]], X_t[obs[:, 1]])
    uv_obs = uv_t + (SIG_IMG * rng.standard_normal(uv_t.shape) if noise else 0)
    gcp_obs = X_t[gcp_ids] + (SIG_GCP * rng.standard_normal((n_gcp, 3)) if noise else 0)
    # start: centres off by 2 m, angles by 1 deg, points by 1 m; nominal intrinsics if self-calibrating
    C0 = C_t + rng.normal(0, 2.0, C_t.shape)
    R0s = Rotation.from_rotvec(np.deg2rad(rng.normal(0, 1.0, (len(R_t), 3)))).as_matrix() @ R_t
    K0 = np.array([2500.0 * 1.02, 1500.0, 1000.0, 0.0, 0.0]) if selfcal else K_TRUE.copy()
    X0 = X_t + rng.normal(0, 1.0, X_t.shape)
    R, C, K, X, cost, its = adjust(R0s, C0, K0, X0, obs, uv_obs, gcp_ids, gcp_obs, selfcal)
    e = X[chk_ids] - X_t[chk_ids]
    eg = X[gcp_ids] - gcp_obs
    return dict(rmse_r=np.sqrt(np.mean(e[:, 0] ** 2 + e[:, 1] ** 2)), rmse_z=np.sqrt(np.mean(e[:, 2] ** 2)),
                mean_z=e[:, 2].mean(), ctrl_3d=np.sqrt(np.mean(np.sum(eg ** 2, 1))),
                dK=K - K_TRUE, cost=cost, its=its, n=n, n_obs=len(obs), ez=e[:, 2], chk=X_t[chk_ids])


GSD = ALT / K_TRUE[0]                                                  # metres per pixel at z = 0
SEEDS = range(5)
R_, C_, tie_, chk_ = build_scene()
o = experiment("4 corners", False, 0, noise=False)
print(f"{len(R_)} nadir images, {o['n']} points ({len(chk_)} check points), {o['n_obs']} image observations,"
      f" GSD {100 * GSD:.1f} cm at z = 0")
print(f"noise-free check (4 corners, fixed): RMSE_r {o['rmse_r'] * 1e3:.1e} mm, RMSE_z {o['rmse_z'] * 1e3:.1e} mm")
for selfcal in (False, True):
    print("\nself-calibration (f, cx, cy, k1, k2 free)" if selfcal else "\ncamera pre-calibrated (intrinsics fixed)")
    print("GCPs                 check RMSE_r cm      check RMSE_z cm      mean dz cm  control RMSE cm"
          + ("  f error px" if selfcal else ""))
    for name in CONFIGS:
        res = [experiment(name, selfcal, s) for s in SEEDS]
        rr = np.array([q["rmse_r"] for q in res]) * 100
        rz = np.array([q["rmse_z"] for q in res]) * 100
        mz = np.mean([q["mean_z"] for q in res]) * 100
        rc = np.mean([q["ctrl_3d"] for q in res]) * 100
        line = (f"{name:18s} {rr.mean():6.2f} [{rr.min():5.2f}-{rr.max():5.2f}]"
                f"  {rz.mean():6.2f} [{rz.min():5.2f}-{rz.max():6.2f}]  {mz:+7.2f}     {rc:6.2f}")
        if selfcal:
            line += f"     {np.mean([q['dK'][0] for q in res]):+6.2f}"
        print(line)

for name in ("4 corners", "4 clustered"):
    q = experiment(name, True, 0)
    print(f"\nvertical error (cm) at check points, self-calibration, {name}, seed 0 (rows: y up)")
    grid = q["ez"].reshape(5, 6) * 100
    for row in grid[::-1]:
        print("  " + " ".join(f"{v:+7.1f}" for v in row))

60 nadir images, 1999 points (30 check points), 10148 image observations, GSD 4.0 cm at z = 0
noise-free check (4 corners, fixed): RMSE_r 2.7e-12 mm, RMSE_z 1.2e-11 mm

camera pre-calibrated (intrinsics fixed)
GCPs                 check RMSE_r cm      check RMSE_z cm      mean dz cm  control RMSE cm
4 corners            1.50 [ 1.27- 1.80]    2.51 [ 1.84-  3.29]    -0.44       0.79
5 corners+centre     1.38 [ 1.29- 1.51]    2.92 [ 1.87-  5.10]    -1.20       0.79
9 spread (3x3)       1.40 [ 1.25- 1.65]    2.38 [ 1.80-  2.96]    +0.87       0.86
16 spread (4x4)      1.18 [ 1.08- 1.32]    2.26 [ 1.87-  2.48]    -0.04       1.06
4 clustered          4.44 [ 2.30- 7.60]    5.39 [ 3.29-  7.37]    -2.35       0.58

self-calibration (f, cx, cy, k1, k2 free)
GCPs                 check RMSE_r cm      check RMSE_z cm      mean dz cm  control RMSE cm  f error px
4 corners            1.81 [ 1.40- 2.40]   13.79 [ 2.29- 32.16]   -11.74       0.65      -1.08
5 corners+centre     1.41 [ 1.32- 1.51]    3.17 [ 1.90-  4.77]    -1.21       0.69      +2.02
9 spread (3x3)       1.42 [ 1.24- 1.63]    2.87 [ 1.80-  4.65]    +0.97       0.82      +1.10
16 spread (4x4)      1.19 [ 1.10- 1.33]    2.41 [ 2.21-  2.50]    -0.14       1.04      +1.10
4 clustered          5.45 [ 4.77- 6.25]   32.46 [ 5.87- 56.50]    +6.62       0.55      +2.57

vertical error (cm) at check points, self-calibration, 4 corners, seed 0 (rows: y up)
     -0.0    -1.1    -2.6    +3.1    -0.9    -0.7
     -0.2    -0.3    +0.3    +3.1    +2.4    -0.2
     -2.6    +1.9    +0.7    +1.6    +4.8    -0.1
     -0.7    +0.6    +5.1    +0.8    +1.1    -3.0
     -5.2    -2.1    +0.9    +1.0    +1.3    -3.5

vertical error (cm) at check points, self-calibration, 4 clustered, seed 0 (rows: y up)
    -48.9   -47.4   -56.9   -75.1   -89.9  -119.2
    -26.7   -25.8   -33.9   -52.0   -71.2  -103.5
     -9.1    -9.5   -20.7   -34.6   -58.4   -93.0
     -0.1    -1.9   -11.8   -26.0   -51.0   -85.0
     +2.8    +2.7    -4.7   -22.5   -42.3   -75.4

A calibrated camera and spread control. With the intrinsics fixed, every configuration whose control surrounds the block gives check-point errors of about a third of the GSD horizontally and about half to three quarters of it vertically, and adding control from four to sixteen points improves the horizontal RMSE and narrows the spread over seeds but changes the averages only modestly; in our reading, the image noise on the tie points, not the control, then limits accuracy. The vertical error is 1.7 to 2.1 times the horizontal, the same direction as the factor of 2.5 Sanz-Ablanedo et al. (2018) report from their real block, though the set-ups are not comparable. Four points clustered in one corner are three times worse horizontally and twice vertically, because the rest of the block is extrapolated from a small control base.

Self-calibration needs the control spread out. With five or more spread points, freeing the camera costs little: the check-point RMSEs are within a few millimetres of the fixed-camera ones. With control only at the four corners, the vertical RMSE averages 13.8 cm, more than three GSD, and ranges from 2.3 to 32 cm over the five seeds, with a mean vertical error of −11.7 cm. In our reading, nothing inside the block holds the camera parameters and the block’s height together, so the noise of each seed decides how far they drift. Seed 0, whose grid is printed, happens to be one of the good cases, with errors of a few centimetres and no visible pattern; the printed output does not show the shape of the error for the other seeds, and our reading that it is the coupling of interior orientation and height that Sanz-Ablanedo et al. describe, rather than a pure offset, is an interpretation. A single control point in the centre removes most of it (3.2 cm). With the four clustered points, the printed grid shows the vertical error growing smoothly with distance from the control corner (bottom left) to about 1.2 m at the far corner; the errors do not change at a constant rate along the rows, so in our reading the surface is bent as well as tilted.

The control points say nothing about it. The RMSE at the control points stays between 0.55 and 1.06 cm in every row, near the 1 cm survey noise, and it is lowest for the clustered configuration, the least accurate one: a few control points are easy for the adjustment to fit, exactly as Sanz-Ablanedo et al. (2018) found in their control-against-check comparison. It rises slightly as control is added, which they also observed and explain by the network becoming more rigid. Only the independent check points, spread over the whole block and more than the NSSDA’s minimum of 20, expose the difference between a 2 cm and a 30 cm survey. The example’s block is a weak one by design, a single nadir grid with no oblique or cross-strip images. James et al. (2017) found self-calibrated and fixed cameras equally accurate by check-point RMSE with as few as four control points in strong distributions on their real gully survey, though their surface models showed systematic differences that the RMSE concealed; the simulated block differs from theirs, and its corner-only result varies widely from seed to seed, so the two are not directly comparable.

Applications

Topographic mapping and change. Aerial photogrammetry produced topographic elevation maps from overlapping images, first with stereo plotters and then with automated matching (Szeliski, 2022, introduction to Chapter 12), and today supplies surface models and orthophotos, which lidar measures actively; satellite stereo is part of remote sensing. Repeat surveys measure change: James et al. (2017) aimed at sub-decimetre change for gully erosion in Morocco and at the velocity of the Super-Sauze landslide in the French Alps, where their processing halved the planimetric error to about 0.1 m and, by their account, doubled the frequency at which velocity changes could be detected. Change detection covers the image-based side.

Engineering geology and mining. Giordan et al. (2020) review drone photogrammetry for landslides, rock walls, debris flows and other engineering-geology applications, and Sanz-Ablanedo et al. (2018) surveyed an open-cast coal mine. Giordan et al. cite a workflow by Nex and Remondino (2014) that allots 15% of the effort to measuring GCPs in the field, and add that with traditional photogrammetry surveying the GCPs usually takes twice or more the time of the drone survey itself, which is why the number of control points is a cost as well as an accuracy question.

Archaeology and cultural heritage. Remondino et al. (2011) list 3D documentation and mapping of archaeological sites and structures among the application domains of low-altitude image-based surveys and describe three projects: the block at Veio, whose orientation used over 330,000 automatically extracted correspondences reduced to about 18,000 for the adjustment and whose MicMac matching produced about 40 million points; the two main courts of the Maya site of Copán in Honduras, where dense vegetation caused matching blunders while non-occluded structures were reconstructed correctly; and an excavation at Pava, Italy, surveyed at the start and end of each excavation season to monitor progress, compute excavated volumes and produce multi-temporal orthoimages. For objects with a prevalent vertical extension, such as a large building façade, and where satellite visibility is limited, Remondino et al. note that direct georeferencing is unsuitable and orientation has to rely on the images alone.

Engineering and industrial measurement. Close-range photogrammetry measures structures, deformations and parts with convergent networks, coded targets and calibrated cameras; the triangulation article describes calibrated multi-camera systems that triangulate targets in every frame, and the direct linear transform article the close-range origin of that method. Single images can also be measured when the scene supplies geometry, as in the single-view metrology of vanishing points.

Computer vision. Photogrammetric reconstructions are also data: SfM pipelines label training sets, and benchmarks such as ETH3D and Tanks and Temples register images to laser scans to obtain independent reference geometry.

Software

The sources read here used both commercial and open-source software; this section describes open-source packages only as their official documentation does.

• COLMAP (its article) implements incremental structure from motion and multi-view stereo; for surveying it offers geo-registration by a RANSAC-estimated similarity from camera-centre coordinates, and mapping with GPS position priors (COLMAP 4.2.1 documentation, as of November 2026).
• OpenDroneMap turns drone images into georeferenced point clouds, meshes, orthophotos and elevation models, as aerial and drone image analysis summarizes; its camera option against doming is described above.
• MicMac is described by its repository as free open-source photogrammetric software for 3D reconstruction, developed at IGN, the French mapping agency, and its national school of geographic sciences (ENSG) within the LASTIG laboratory, and distributed under the CeCILL-B licence since 2007; a second version, MMVII, is developed in a separate repository (README, as of November 2026). Remondino et al. (2011) used its matching for the Veio surface.
• The geoscience studies above used Agisoft PhotoScan (James et al., 2017, version 1.1.6; Sanz-Ablanedo et al., 2018). James et al. (2017) note that PhotoScan gives no algorithmic details of its “optimization” and does not report the precisions and covariances of the camera parameters, which is why they built the Monte Carlo approach around it.

Challenges

• Weak network geometry. An accurate camera model matters most with the relatively parallel and flat image geometries typical of aerial surveys (James et al., 2017), and correlation between inner and exterior orientation is the condition behind doming (Sanz-Ablanedo et al., 2018); James et al. (2017) also note that weak or noisy networks, particularly with strongly distorting lenses, need more control.
• Unstable consumer cameras. Lens settings can change between flights, as in the landslide survey of James et al. (2017), where EXIF data showed two focal lengths and two camera models had to be used.
• Ground control that is wrong or moves. Misidentified or disturbed targets and moving ground bias the adjustment, as the GCP removed for tilting the model and the 47 removed on the landslide show (James et al., 2017).
• Targets that are hard to see. At 6.86 cm GSD the 80 cm targets of Sanz-Ablanedo et al. (2018) were about 10 pixels across, and those on low ground appeared smaller and were harder to measure; Giordan et al. (2020) warn that white markers can saturate.
• Surfaces that do not match. Vegetation produced matching blunders at Copán (Remondino et al., 2011) and was likely partly responsible for large tie-point residuals on the landslide (James et al., 2017); snow, sand and water are discussed in aerial and drone image analysis, and weak texture in multi-view stereo.
• Cost of control. Surveying GCPs is a large share of field effort (Giordan et al., 2020), which is the practical pressure towards fewer points and direct georeferencing.

Failure Modes and Risks

• Undetected doming. A dome far from the control points does not raise the control RMSE, and check points placed only where the dome is small do not reveal it either; Sanz-Ablanedo et al. (2018) therefore ask for many check points distributed over the whole area, not only at its edge.
• Accuracy claimed from control residuals. Reporting the fit at control points as accuracy overstates it, by large factors when control is sparse (Sanz-Ablanedo et al., 2018).
• Wrong weights. Overweighted control bends the model towards errors in the control; underweighted control fails to correct systematic error (James et al., 2017).
• Extrapolation beyond the control. Accuracy measured inside a controlled area does not carry to the area outside it (James et al., 2017; the rolling-shutter study described in aerial and drone image analysis found the same).
• Change that is an artefact. Systematic errors that differ between two epochs appear as surface change; on the gully survey of James et al. (2017), two surface models from the default processing, with and without the check points, differed by volumes of 173 m³ on the western side and 676 m³, of opposite sign, on the eastern side.

Deployment Considerations

Reporting. James et al. (2019) give a checklist for publishing SfM photogrammetry: the camera (model, sensor and image size, focal length, lens type and, for moving platforms, shutter type); the survey design and its rationale; the number, distribution, measurement method and assumed precision of the control, with check points identified separately; deviations from the design and field conditions; the software with its version and all relevant processing settings, including which camera parameters were self-calibrated and which control entered the adjustment; image residuals at tie and control points; accuracy from independent check data, with bias, precision and spatial variability; and the propagation of uncertainty into derived products.

Standards. Positional accuracy is reported against a stated standard, such as the NSSDA statistics above, and the thresholds are set by the agency, contract or application, not by the standard (Federal Geographic Data Committee, 1998). Rules on flying aircraft and drones, which vary by country, are outlined from official texts in aerial and drone image analysis.