Building an Image Pair Graph from GPS Positions

Feature matching over every pair of a large survey is not slow, it is infeasible — fifty million comparisons for ten thousand images. The single most effective reduction uses information the survey already carries: two frames taken far apart cannot overlap, so their pair never needs testing.

Doing that correctly takes more than a distance threshold. The positions are in degrees and the threshold is in metres, the right radius depends on the altitude and the sensor, terrain relief changes the footprint, and a radius that is slightly too small fragments the reconstruction in a way that looks like a matching failure. This page covers each of those, as the implementation detail behind matching strategies for large and linear surveys.

Working in a metric frame

EXIF positions are latitude and longitude in degrees, and a degree is not a distance. A degree of longitude is 111 km at the equator and 71 km at 50° north, so a threshold in degrees is a different threshold at every latitude and a different one in each axis.

Projecting to a local metric frame before doing any geometry removes the problem entirely and costs one transformation. An east-north-up frame centred on the survey is the natural choice: distances are metres, the axes are orthogonal, and nothing about the arithmetic depends on where the site is.

import numpy as np
from pyproj import CRS, Transformer


def to_local_enu(lat: np.ndarray, lon: np.ndarray, alt: np.ndarray) -> dict:
    """Project geographic positions into a local metric frame.

    A transverse Mercator centred on the survey keeps distortion under a part
    per million across any site a drone can fly, so distances in the frame
    are metres to well within the accuracy anything downstream needs.
    """
    lat0, lon0 = float(np.mean(lat)), float(np.mean(lon))
    local = CRS.from_proj4(
        f"+proj=tmerc +lat_0={lat0} +lon_0={lon0} +k=1 +x_0=0 +y_0=0 "
        "+ellps=WGS84 +units=m +no_defs")
    transformer = Transformer.from_crs(CRS.from_epsg(4326), local, always_xy=True)
    east, north = transformer.transform(lon, lat)
    return {"enu": np.column_stack([east, north, alt]),
            "origin": (lat0, lon0), "crs": local.to_proj4()}

Deriving the radius from the geometry

The radius should be the distance beyond which two frames cannot share ground, which follows from the sensor and the altitude rather than from experience.

import numpy as np


def footprint_width(altitude_m: float, sensor_width_mm: float,
                    focal_mm: float) -> float:
    """Ground width covered by one frame, from the sensor geometry."""
    return altitude_m * sensor_width_mm / max(focal_mm, 1e-9)


def pair_radius(altitude_m: float, sensor_width_mm: float, focal_mm: float,
                *, relief_m: float = 0.0, margin: float = 1.15) -> dict:
    """Radius within which two frames may overlap.

    Relief matters in the direction that surprises people: ground closer to
    the aircraft has a smaller footprint, so a hilltop reduces the overlap
    and the radius must be computed at the lowest altitude above ground
    anywhere in the survey rather than at the nominal one.
    """
    effective_altitude = max(altitude_m - relief_m, 1.0)
    width = footprint_width(effective_altitude, sensor_width_mm, focal_mm)
    return {"footprint_m": width, "radius_m": width * margin,
            "effective_altitude_m": effective_altitude,
            "note": ("relief reduces the footprint on high ground, so the radius is "
                     "computed there" if relief_m > 0 else "flat-terrain assumption")}

Computing at the lowest height above ground is the conservative choice and the right one. A survey flown at 100 m over a site with a 40 m hill has a footprint on the hilltop that is only 60 % of the one in the valley, and a radius derived from the nominal altitude will exclude genuine overlaps exactly where the terrain is most interesting.

How terrain relief shrinks the frame footprint and the overlap A cross-section showing an aircraft flying a constant barometric altitude over terrain with a forty metre hill. Over the valley the height above ground is one hundred metres and the frame footprint is one hundred and twenty metres wide, giving generous overlap between consecutive frames. Over the hilltop the height above ground is sixty metres and the footprint is seventy-two metres, so consecutive frames overlap far less. A note states that a pair radius derived from the nominal altitude excludes genuine overlaps exactly where the terrain is most complex. constant flight altitude terrain 120 m footprint 72 m footprint 40 m hill The footprint shrinks where the ground rises, and so does the overlap. A radius from the nominal altitude excludes pairs exactly where the terrain is most complex.

Figure 1 — Why the radius is computed at the lowest height above ground rather than the nominal altitude.

Building the graph

import numpy as np
from scipy.spatial import cKDTree


def build_pair_graph(enu: np.ndarray, radius_m: float,
                     *, sequential_window: int = 4,
                     max_pairs_per_image: int = 400) -> dict:
    """Candidate pairs from position, plus capture-order neighbours.

    Capping per image rather than globally keeps the graph even: a global cap
    truncated in index order would leave the end of a survey unmatched, which
    is a failure mode that produces a reconstruction missing its last lines.
    """
    tree = cKDTree(enu[:, :2])
    pairs = set()

    for i in range(len(enu)):
        neighbours = tree.query_ball_point(enu[i, :2], radius_m)
        neighbours = [j for j in neighbours if j != i]
        if len(neighbours) > max_pairs_per_image:
            distances = np.linalg.norm(enu[neighbours, :2] - enu[i, :2], axis=1)
            keep = np.argsort(distances)[:max_pairs_per_image]
            neighbours = [neighbours[k] for k in keep]
        for j in neighbours:
            pairs.add((min(i, j), max(i, j)))

    for i in range(len(enu)):
        for j in range(i + 1, min(i + sequential_window + 1, len(enu))):
            pairs.add((i, j))

    return {"pairs": sorted(pairs), "count": len(pairs),
            "mean_per_image": 2 * len(pairs) / max(len(enu), 1)}
What a position-derived pair graph must include beyond nearest neighbours Three rows. Along-track neighbours are the obvious pairs and are what a naive radius query returns, but a graph containing only these produces a chain that drifts, because nothing ties one flight line to the next. Cross-track pairs between adjacent lines are what make the block rigid, and they are the pairs a simple distance threshold under-weights because side overlap is smaller than forward overlap. Loop closures, where the flight revisits an earlier area, are few in number and disproportionately valuable, because they are the only observations that constrain accumulated drift. along-track neighbours what a naive radius query returns — a chain that drifts cross-track pairs what makes the block rigid; under-weighted by a plain distance threshold loop closures few, and the only observations that constrain accumulated drift A graph of along-track pairs alone reconstructs beautifully and is bent.

Figure 3 — Three pair types, and only the first is found by accident.

Edge-case matrix

Situation Effect Handling
Threshold in degrees Different radius per latitude and axis Project to a metric frame
Nominal altitude used with relief Pairs missing on high ground Compute at the lowest height above ground
Missing positions on some frames Those frames get no spatial pairs Sequential window covers them
Wrong positions Wrong pairs both ways Sanity-check positions first
Very dense capture Hundreds of neighbours per image Cap per image, nearest first
Two flights merged Spurious long-distance pairs Split by time gap before pairing
Variable altitude survey One radius is wrong everywhere Per-image radius from each frame’s altitude
Oblique orbit around a structure Positions close, views opposed Add appearance candidates

The oblique-orbit row is a genuine limitation of position-based pairing. Two frames on opposite sides of a building are metres apart and share no ground at all, while two frames on the same side are further apart and overlap completely. Position alone cannot tell them apart, and the remedy is appearance-based candidates or gimbal-attitude filtering.

import numpy as np


def filter_by_view_direction(enu: np.ndarray, yaw_deg: np.ndarray,
                             pairs: list, *, max_angle_deg: float = 100.0) -> list:
    """Drop pairs whose cameras were pointing in incompatible directions.

    Two frames a few metres apart looking in opposite directions share no
    ground, and testing the pair wastes matching time and occasionally
    produces a confident wrong match on a repetitive facade.
    """
    yaw = np.radians(np.asarray(yaw_deg, dtype=float))
    direction = np.column_stack([np.cos(yaw), np.sin(yaw)])
    limit = np.cos(np.radians(max_angle_deg))
    return [(i, j) for i, j in pairs
            if float(direction[i] @ direction[j]) > limit]

Verification snippet

import networkx as nx
import numpy as np


def verify_pair_graph(pairs: list, n_images: int,
                      *, min_median_degree: int = 8) -> dict:
    """Connectivity and redundancy checks before any matching runs."""
    graph = nx.Graph()
    graph.add_nodes_from(range(n_images))
    graph.add_edges_from(pairs)

    components = list(nx.connected_components(graph))
    degrees = np.array([d for _, d in graph.degree()])
    problems = []

    if len(components) > 1:
        sizes = sorted((len(c) for c in components), reverse=True)
        problems.append(f"{len(components)} components, largest {sizes[0]} of {n_images} "
                        "— the reconstruction will fragment")
    if float(np.median(degrees)) < min_median_degree:
        problems.append(f"median degree {np.median(degrees):.0f} is below "
                        f"{min_median_degree} — single match failures will disconnect images")
    isolated = int((degrees == 0).sum())
    if isolated:
        problems.append(f"{isolated} images have no candidate pairs at all")

    return {"components": len(components), "median_degree": float(np.median(degrees)),
            "min_degree": int(degrees.min()), "isolated": isolated,
            "problems": problems, "ready": not problems}

Running this before matching converts hours of wasted compute into a message. A disconnected graph will always produce a fragmented reconstruction, and it takes a second to detect and minutes to fix by widening the radius.

Median graph degree and reconstruction outcome against pair radius Two curves against the pair radius expressed as a multiple of the frame footprint, from zero point six to two times. Median graph degree rises from three at zero point six to twenty-six at two times. Reconstruction success, shown as the fraction of images included, is only sixty percent at zero point six, rises steeply to ninety-nine percent by one point one and is flat thereafter. A knee is marked at one point one five, where connectivity is achieved and the pair count has not yet begun to grow quadratically for nothing. 0.6×0.9× 1.15×1.5×2.0× pair radius, as a multiple of the frame footprint median degree images reconstructed the knee — 99 % at the smallest radius Beyond the knee the pair count keeps growing quadratically and the outcome does not improve.

Figure 2 — Where to choose the radius, measured on the survey rather than assumed.

Splitting merged flights before pairing

A directory containing two flights produces spatial pairs between them wherever they overlap, which is usually correct, and sequential pairs across the join, which is not: the last frame of one flight and the first of the next are adjacent in index order and may be kilometres apart.

import numpy as np


def split_by_time_gap(timestamps_s: np.ndarray, *, gap_s: float = 600.0) -> list[dict]:
    """Group frames into flights by looking for long gaps in capture time."""
    order = np.argsort(timestamps_s)
    t = np.asarray(timestamps_s, dtype=float)[order]
    breaks = np.flatnonzero(np.diff(t) > gap_s) + 1
    groups = np.split(order, breaks)
    return [{"flight": i, "indices": g.tolist(),
             "duration_min": float((t[breaks[i - 1] if i else 0] - t[0]) / 60)}
            for i, g in enumerate(groups)]

Sequential pairing should then run within each flight rather than across the whole directory, which is a one-line change and removes a small population of nonsensical pairs that occasionally produce a confident wrong match.

When to escalate

  • The graph is disconnected at every radius. Coverage has a genuine gap — the survey missed an area, or a section was flown at a different altitude. Widening will not bridge it.
  • Positions are present and the graph still fragments. Check the positions first; a hemisphere or datum error produces exactly this.
  • The survey has no positions. Fall back to appearance-based pairing, covered in using vocabulary trees for large dataset matching.

← Matching Strategies for Large and Linear Surveys