Orthomosaic Radiometry and Seamline Control

The checkpoints pass at two centimetres. The DEM is clean. And the client’s first email says the orthomosaic looks like a patchwork quilt — visible rectangles where flight lines meet, a bright stripe down the middle of the site, half a lorry on a haul road.

None of that is a geometry problem, which is why none of it shows up in an accuracy report. It is radiometry: the brightness, colour and consistency of the pixels, and the choice of where in the overlap the pipeline cut from one image to the next. This page covers that half of orthomosaic quality — why the brightness varies at all, what a seamline is doing, how to make the cut fall somewhere the eye will not follow, and how to measure the result rather than squint at it.

It sits within DEM/DSM generation and raster export automation, downstream of reconstruction and upstream of exporting Cloud-Optimized GeoTIFFs with rasterio.

Audience and prerequisites. Python 3.10+, rasterio and NumPy, an orthomosaic you can regenerate, and access to the source images with their EXIF exposure fields intact.

# Prerequisites

# Conceptual architecture

Four independent causes push two overlapping images to disagree about the colour of the same patch of ground. Treating them as one problem — “the mosaic looks wrong” — is why the usual fix, a global contrast stretch at the end, never works.

Exposure variation. An auto-exposing camera changes shutter and ISO between frames. A flight line flown into brightening sky produces a ramp of exposure across a dozen images.

Illumination change. The sun moves during a flight, and cloud moves faster. A forty-minute survey can have two quite different lighting conditions in it, and the boundary between them will sit somewhere in the middle of the site.

Bidirectional reflectance. The ground genuinely reflects different amounts of light in different directions. The half of each image facing away from the sun is darker than the half facing it — this is the hotspot, and it survives every exposure correction because it is real.

Vignetting and lens falloff. Every frame is darker at its corners than its centre, by a fixed amount determined by the optics.

Figure 1 — Four causes, four scales, four different corrections.

The order matters. Vignetting and exposure are per-image and are corrected before anything is compared. Only then does a block adjustment of residual gains make sense, because otherwise it is trying to absorb a lens effect into a per-image constant and cannot.

# The seamline, and why it is a choice

Where two corrected images still disagree slightly, the mosaic has to pick one. The line along which it switches is the seamline, and the mosaicking step has enormous freedom in where to put it.

A naive mosaic cuts on the footprint boundary — a straight line through the middle of the overlap. That is the worst possible choice, because a straight cut across a uniform field makes any residual difference into a visible edge that runs for hundreds of metres. The eye is extremely good at detecting a straight intensity step and almost blind to the same step following a hedge.

A good seamline does three things:

• Follows existing edges. A cut along a field boundary, a kerb, a fence or a tree line hides a residual difference inside a discontinuity the viewer already expects.
• Avoids tall objects. A seamline crossing a building or a mast produces a sheared roof, because the two images see the tall object leaning in different directions.
• Stays away from moving things. A cut through a car park will duplicate or truncate vehicles.

Figure 2 — The same residual difference, placed well and placed badly.

# Minimal reproducible solution

Before changing anything, measure. This function gives the two numbers that tell you whether the problem is exposure, seamlines or neither.

import numpy as np
import rasterio
from rasterio.windows import Window


def mosaic_radiometry_report(path: str, *, block: int = 512) -> dict:
    """Per-block brightness statistics across a mosaic.

    Two signals matter. A wide spread of block means says the images were not
    equalised before mosaicking. A narrow spread with a high maximum gradient
    between neighbouring blocks says the equalisation worked and the seamlines
    are where the remaining visible problem is.
    """
    means, positions = [], []
    with rasterio.open(path) as src:
        nodata = src.nodata
        for row in range(0, src.height, block):
            for col in range(0, src.width, block):
                win = Window(col, row,
                             min(block, src.width - col),
                             min(block, src.height - row))
                data = src.read(window=win, masked=True).astype("float32")
                if data.count() < 0.5 * data.size:
                    continue                      # mostly nodata; skip the edge
                means.append(float(data.mean()))
                positions.append((row // block, col // block))

    arr = np.array(means)
    grid = {p: m for p, m in zip(positions, arr)}
    gradients = [
        abs(grid[(r, c)] - grid[(r, c + 1)])
        for (r, c) in grid if (r, c + 1) in grid
    ] + [
        abs(grid[(r, c)] - grid[(r + 1, c)])
        for (r, c) in grid if (r + 1, c) in grid
    ]

    return {
        "blocks": len(arr),
        "mean": float(arr.mean()),
        "spread_pct": float(100 * arr.std() / arr.mean()),
        "p95_gradient": float(np.percentile(gradients, 95)) if gradients else 0.0,
        "max_gradient": float(max(gradients)) if gradients else 0.0,
        "nodata": nodata,
    }

A spread above about 8% means the per-image equalisation did not happen or did not work. A spread under 4% with a 95th-percentile block-to-block gradient above roughly 6 digital numbers means equalisation is fine and seamlines are the issue. Both high means both.

# Parameter deep-dive

Blending width. Feathering across a transition zone hides a step by spreading it. Too narrow and the step is merely blurred; too wide and any misregistration turns into a double image. A width of 30–80 pixels works for most survey mosaics. It is a cosmetic fix and it does not repair a large radiometric difference — it converts a sharp seam into a soft gradient, which is better but still visible against uniform ground.

Gain versus offset. Correcting brightness by a multiplicative gain preserves contrast ratios and is right for exposure differences. An additive offset is right for haze and path radiance. Using the wrong one is a common cause of a mosaic that matches in the midtones and diverges in the shadows.

Reference image selection. A block adjustment needs an anchor. Choosing the first image alphabetically means the whole mosaic inherits the exposure of one arbitrary frame. Choosing the image closest to the block’s median brightness leaves the smallest total correction and the least clipping.

Overlap sampling. Gains are estimated from pixels seen by both images. Sampling the whole overlap includes tall objects and moving vehicles, which bias the estimate. Sampling only the flat, static parts — or trimming the top and bottom deciles of the difference distribution — is markedly more robust.

Clipping policy. A gain above about 1.3 applied to an 8-bit image will clip highlights permanently. Capping gains and accepting a small residual mismatch is nearly always better than destroying detail in bright areas.

# Estimating per-image gains as a block adjustment

Correcting each image against its neighbour in isolation does not work: the corrections chain, and a long flight line accumulates a drift that ends with the last image several stops away from the first. The problem is global, and it has the same shape as a levelling network — a set of observed differences between pairs, and one unknown per image, solved together.

Each overlapping pair contributes one observation: the ratio of mean brightness over the pixels both images see. Taking logarithms turns the multiplicative gain into an additive unknown, which makes the whole thing an ordinary least-squares problem.

import numpy as np


def solve_block_gains(pairs: list[tuple[int, int, float]], n_images: int,
                      *, anchor: int = 0) -> np.ndarray:
    """Solve for one gain per image from pairwise brightness ratios.

    `pairs` holds (image_a, image_b, ratio) where ratio is mean(a) / mean(b)
    over the pixels the two images share. Working in log space turns each
    observation into g_a - g_b = log(ratio), a linear equation, and the system
    is rank-deficient by exactly one — the overall brightness of the mosaic is
    unconstrained — so one image is anchored.
    """
    rows, obs = [], []
    for a, b, ratio in pairs:
        if ratio <= 0:
            continue
        row = np.zeros(n_images)
        row[a], row[b] = 1.0, -1.0
        rows.append(row)
        obs.append(np.log(ratio))

    # The anchor equation: g_anchor = 0, i.e. that image keeps its own exposure.
    row = np.zeros(n_images)
    row[anchor] = 1.0
    rows.append(row)
    obs.append(0.0)

    A = np.vstack(rows)
    log_gains, *_ = np.linalg.lstsq(A, np.array(obs), rcond=None)
    return np.exp(log_gains)


def choose_anchor(mean_brightness: np.ndarray) -> int:
    """Anchor on the image nearest the block median, not the first one.

    Anchoring arbitrarily makes every other image move toward one frame's
    exposure; anchoring on the median leaves the smallest total correction and
    the least clipping at both ends of the range.
    """
    return int(np.argmin(np.abs(mean_brightness - np.median(mean_brightness))))

Two practical notes. Pairs with very small overlap should be dropped — a ratio computed over two hundred pixels is noise, and it will drag the solution. And the residuals of this fit are diagnostic in their own right: a pair whose observed ratio disagrees strongly with the solved gains is usually a pair where one image caught a cloud shadow, which is worth knowing before the mosaic is built rather than after.

# Routing a seamline through the overlap

Once gains are solved, the remaining difference has to be put somewhere. Treating the overlap as a grid where every pixel has a cost, and finding the cheapest path across it, is the standard formulation — and the cost function is where all the judgement lives.

Three terms make up a cost that behaves well in practice:

Difference cost. The absolute difference between the two images at that pixel. Cutting where they already agree is free; cutting where they disagree by twenty digital numbers is expensive.

Edge bonus. The local gradient magnitude of the image. A strong edge reduces the cost, because a cut hidden in a hedge or a kerb is invisible. This is the term that produces the behaviour in Figure 2.

Exclusion cost. A very large constant over tall objects and moving vehicles, which makes the path route around them rather than through.

from scipy import ndimage


def seam_cost_surface(img_a: np.ndarray, img_b: np.ndarray,
                      exclusion: np.ndarray, *,
                      edge_weight: float = 0.6,
                      exclusion_cost: float = 1e6) -> np.ndarray:
    """Build the cost grid a seamline search will traverse.

    The edge term is subtracted rather than added: a strong image gradient is
    somewhere a cut can hide, so it should attract the path, not repel it.
    """
    difference = np.abs(img_a.astype("float32") - img_b.astype("float32"))

    blended = 0.5 * (img_a.astype("float32") + img_b.astype("float32"))
    gx = ndimage.sobel(blended, axis=1)
    gy = ndimage.sobel(blended, axis=0)
    edges = np.hypot(gx, gy)
    edges = edges / (edges.max() or 1.0)

    cost = difference - edge_weight * difference.max() * edges
    cost = np.clip(cost, 0.0, None) + 1.0          # keep every step positive
    cost[exclusion.astype(bool)] = exclusion_cost
    return cost

The exclusion mask is worth building properly. A height raster thresholded a metre or two above the local ground plane catches buildings, masts and mature trees; a difference mask between the two images, thresholded high and then dilated, catches most moving vehicles. Both are cheap to compute and between them they prevent the two most conspicuous mosaic defects.

# Edge-case matrix

The last row causes more confusion than it should. A block statistic that silently includes the black border of the mosaic produces a spread figure that looks alarming and means nothing.

# Verification snippet

Measuring across the seams specifically is what separates “it looks fine to me” from a check you can put in a pipeline.

from scipy import ndimage


def seam_step_magnitude(mosaic_path: str, seam_mask_path: str,
                        *, band: int = 1, buffer_px: int = 4) -> dict:
    """Measure the intensity step across seamlines rather than everywhere.

    A seam mask is a raster of the same shape marking seamline pixels. Dilating
    it by a few pixels on each side and comparing the two sides gives the step
    magnitude that a viewer actually sees, which a whole-image statistic dilutes
    into invisibility.
    """
    with rasterio.open(mosaic_path) as src:
        img = src.read(band, masked=True).astype("float32")
    with rasterio.open(seam_mask_path) as src:
        seam = src.read(1) > 0

    near = ndimage.binary_dilation(seam, iterations=buffer_px)
    ring = near & ~seam
    labels, n = ndimage.label(ring)
    if n < 2:
        return {"seams": 0, "mean_step": 0.0, "max_step": 0.0}

    steps = []
    for start in range(1, n, 2):
        a = img[labels == start].mean()
        b = img[labels == start + 1].mean() if start + 1 <= n else a
        steps.append(abs(float(a) - float(b)))

    return {
        "seams": len(steps),
        "mean_step": float(np.mean(steps)),
        "max_step": float(np.max(steps)),
        # Below about 3 DN on 8-bit imagery the seam is effectively invisible
        # against anything but a perfectly uniform surface.
        "acceptable": float(np.max(steps)) < 3.0,
    }

Figure 3 — Order of operations. Each stage assumes the previous one has run.

# FAQ

# Can a bad mosaic be fixed without re-running the whole reconstruction?

Often, yes. Radiometry and seamlines are applied at the mosaicking stage, which is the last step and usually the cheapest to repeat. If the reconstruction is sound, re-mosaicking with corrected inputs and better seamline routing costs a fraction of a full re-run — which is a good reason to keep the aligned project rather than only the exported deliverables.

# Why does the mosaic look fine at full extent and terrible when zoomed in?

Overviews average away small steps. A mosaic viewed at 1:5000 through a decimated pyramid level can look seamless while the full-resolution pixels show a clear edge. Always inspect at native resolution, and treat a check that only runs on overviews as no check at all.

# Does this matter for a DEM as well as an orthomosaic?

Radiometry does not, because a DEM has no colour. Seamlines very much do: a height discontinuity where two surfaces meet is the elevation equivalent of a visible seam, and it propagates into contours and volumes. The surface generation work handles that side.

# Should multispectral mosaics use the same approach?

No. Equalising gains to make a mosaic look uniform destroys exactly the reflectance differences a multispectral survey exists to measure. Those mosaics are calibrated to reflectance per band instead — see radiometric calibration in Python.

# Guides in this topic

• Colour Balancing Images Before Mosaicking
• Fixing Visible Seamlines and Colour Banding
• Choosing Blending Modes for Orthomosaic Quality
• Removing Moving Objects and Ghosting from Mosaics

# Related

• Exporting Cloud-Optimized GeoTIFFs with Rasterio
• Publishing Rasters to STAC and Object Storage
• Troubleshooting DEM and Raster Export Failures

← DEM/DSM Generation & Raster Export Automation