Automating Camera Intrinsic Matrix Extraction

You run a batch of imagery from a mixed UAV fleet through Structure-from-Motion and the reconstruction either refuses to converge or produces an orthomosaic that is visibly stretched, sheared, or scaled wrong — yet the same pipeline worked perfectly last week on a single airframe. The usual culprit is the camera intrinsic matrix (K): a focal length transcribed by hand, a value scraped from an unverified EXIF dump, or a single hardcoded K reused across two payloads with different sensors. When K is wrong, every 2D-to-3D projection is wrong, so bundle adjustment diverges and georeferencing drifts. This page shows how to extract and validate K automatically at ingest, before a single feature is detected, so a swapped gimbal or a firmware revision can no longer silently corrupt a run.

# Why intrinsic errors surface downstream, not at ingest

The intrinsic matrix defines the optical projection geometry of the payload — it maps 3D scene rays onto the 2D pixel grid. UAV manufacturers do not store K directly; they store focal length in millimetres and sensor dimensions, and the pipeline is expected to derive pixel-space focal lengths (fx, fy) itself. That derivation is deterministic, but the inputs are fragile, which is why a bad K rarely throws an error at the point it is built and instead detonates later inside the solver.

The conversion from physical millimetres to pixels is:

$$f_x = \frac{f_{\text{mm}}}{w_{\text{sensor}}}\, W_{\text{px}} \qquad f_y = \frac{f_{\text{mm}}}{h_{\text{sensor}}}\, H_{\text{px}}$$

The principal point $(c_x, c_y)$ defaults to the image centre $(W_{\text{px}}/2,\ H_{\text{px}}/2)$ unless a measured calibration value is available. Lens distortion coefficients (k1, k2, p1, p2, k3) live separately in EXIF MakerNotes or manufacturer calibration files (.cal, .xml); they are applied during undistortion and must never be folded into the 3×3 intrinsic matrix.

Three input failures account for nearly every corrupted K in a fleet pipeline:

• Stripped MakerNotes. Aggressive post-processing (DJI Fly exports, Lightroom round-trips, automated cloud uploaders) drops proprietary calibration tags, so the principal point and distortion silently fall back to defaults.
• 35 mm-equivalent focal length misuse. FocalLengthIn35mmFilm is populated for marketing and is wrong for metric calibration — using it inflates fx/fy by the crop factor and guarantees solver divergence.
• Non-square pixel artefacts. Older payloads or transcoded video frames can present fx ≠ fy. Modern UAV cameras hold |fx - fy| / max(fx, fy) < 0.02; a larger spread signals corrupted EXIF or a wrong sensor-dimension override.

Because K is consumed by every later stage, it has to be validated at the boundary. This is the calibration dependency that feeds both the feature detection algorithms for drone imagery that match keypoints across frames and the spatial-reference enforcement covered in managing coordinate reference systems in GDAL, where the projected point cloud must land in the correct datum.

# Minimal reproducible solution

The function below is the core of the extractor: it reads the EXIF focal length and image dimensions, derives fx/fy/cx/cy, applies the validation thresholds, and returns an OpenCV-compatible K plus a structured report. Keep it under one screen; the production wrapper (CLI flags, directory walking, fallback loading) is layered around this without changing the maths.

import os
from typing import Dict, Optional, Tuple

import numpy as np
import exifread
from PIL import Image


def extract_intrinsic_matrix(
    image_path: str,
    sensor_width_mm: float,
    sensor_height_mm: float,
    max_focal_dev_pct: float = 15.0,   # focal symmetry tolerance
    pp_tol_frac: float = 0.05,         # principal-point drift tolerance
    principal_point: Optional[Tuple[float, float]] = None,
) -> Tuple[np.ndarray, Dict]:
    """Derive and validate the 3x3 camera intrinsic matrix from EXIF."""
    with open(image_path, "rb") as f:
        tags = exifread.process_file(f, details=False)

    width_px, height_px = Image.open(image_path).size

    focal_tag = tags.get("EXIF FocalLength")          # real focal, NOT 35mm-equiv
    if focal_tag is None:
        raise ValueError(f"Missing EXIF FocalLength in {os.path.basename(image_path)}")
    focal_mm = float(focal_tag.values[0].num) / float(focal_tag.values[0].den)

    fx = (focal_mm / sensor_width_mm) * width_px      # pixel focal lengths
    fy = (focal_mm / sensor_height_mm) * height_px

    # Use a measured principal point when calibration supplies one, else centre.
    cx, cy = principal_point if principal_point else (width_px / 2.0, height_px / 2.0)

    focal_dev = abs(fx - fy) / max(fx, fy) * 100.0    # validation metrics
    pp_dev = max(abs(cx - width_px / 2.0) / width_px,
                 abs(cy - height_px / 2.0) / height_px)

    report = {
        "image": os.path.basename(image_path),
        "fx": fx, "fy": fy, "cx": cx, "cy": cy,
        "focal_deviation_pct": focal_dev,
        "pp_drift_frac": pp_dev,
        "passed": focal_dev <= max_focal_dev_pct and pp_dev <= pp_tol_frac,
    }

    K = np.array([[fx, 0.0, cx],
                  [0.0, fy, cy],
                  [0.0, 0.0, 1.0]], dtype=np.float64)
    return K, report

The validation thresholds baked into the defaults are the ones that matter in production UAV mapping:

• Focal symmetry: |fx - fy| / max(fx, fy) ≤ 0.15. Exceeding this means non-square pixels, EXIF corruption, or a wrong sensor-dimension override.
• Principal-point drift: |cx - W/2| / W ≤ 0.05 and |cy - H/2| / H ≤ 0.05. Larger offsets indicate severe lens decentering or gimbal misalignment.
• Focal plausibility: fx should fall between 0.8 × W and 1.5 × W. Values outside this band almost always trace back to FocalLengthIn35mmFilm leaking into the calculation.

For automation, wrap this function in an argparse CLI so thresholds and fallbacks are version-controlled rather than edited inline:

# Edge-case matrix

When validation fails, route to fallback strategies in priority order: first a manufacturer calibration file (.cal/.xml bundled with the flight log), then a self-calibration seed (fx = fy = 0.8 × max(W, H), cx = W/2, cy = H/2) handed to cv2.calibrateCamera or a COLMAP PINHOLE model with --ImageReader.single_camera 1, and finally a version-controlled fleet registry mapping CameraModel → K for known payloads.

# Verify the fix worked

After extraction, assert the matrix is structurally sound and physically plausible before letting it reach the solver. This block fails loudly rather than passing a quietly-wrong K downstream:

import numpy as np

K, report = extract_intrinsic_matrix("input/DJI_0001.JPG", 13.2, 8.8)
W = 5472  # native sensor width in px for this payload

assert K.shape == (3, 3) and K[2, 2] == 1.0, "K is not a normalized 3x3"
assert K[0, 1] == 0.0 and K[1, 0] == 0.0, "skew must be zero for UAV pinhole"
assert 0.8 * W <= K[0, 0] <= 1.5 * W, f"fx={K[0,0]:.1f} implausible vs width {W}"
assert report["passed"], f"validation failed: {report}"
print(f"OK fx={K[0,0]:.1f} fy={K[1,1]:.1f} drift={report['pp_drift_frac']:.3%}")

A passing run prints fx/fy within a few percent of each other and a sub-5% principal-point drift; group the per-image reports by EXIF Model and confirm each model resolves to a single stable K before submitting the batch.

# When to escalate

• Reprojection error stays high with a valid K. If bundle adjustment still diverges after K passes every assert, the problem is lens distortion or pose, not intrinsics — move to camera self-calibration and the solver tuning in optimizing bundle adjustment with Python.
• The orthomosaic is correctly shaped but mislocated. A correct K paired with a wrong datum or axis order still drifts on the ground; that is a spatial-reference fault — return to managing coordinate reference systems in GDAL and the survey-grid alignment in coordinate transformation workflows in pyproj.
• Survey-grade absolute accuracy is required. EXIF-derived intrinsics get SfM to converge but are not a substitute for a metric calibration target; tie the result to control as described in setting accuracy thresholds for survey projects.

# Related

• Managing Coordinate Reference Systems in GDAL
• Feature Detection Algorithms for Drone Imagery
• Optimizing Bundle Adjustment with Python
• Structuring Drone Imagery for Batch Processing

← Managing Coordinate Reference Systems in GDAL