Merging Overlapping Orthophoto Tiles Without Seams

Two flights covering the same ground produce two correct images that meet at a visible line, and removing that line is a decision about where to cut rather than a filter to apply afterwards. This page is a focused companion to merge strategies for raster and imagery mosaics.

Concept & Context

When a merge arbitrates a chunk to a different capture than its neighbour, the boundary between them becomes visible. The cause is not registration error — the pixels are in the right place — but radiometry: sun angle, atmospheric scattering and sensor calibration differ between flights, so the same roof is a slightly different colour in each.

Three things reduce the artefact, in increasing order of intrusiveness. Choosing where to cut costs nothing and helps most: a boundary that follows a hedgerow or a river is far less visible than one that crosses a uniform field. Histogram matching brings the two captures into approximate radiometric agreement before anything is blended. Feathering hides what remains, over a band narrow enough not to smear real edges.

The order matters. Feathering a large radiometric difference produces a soft gradient that reads as a smudge rather than a seam, which is not obviously better. Match first, then feather across a narrow band, and the boundary genuinely disappears.

Where the cut line goes decides how visible the seam is Two panes over the same overlap between two captures. In the first, a straight cut along the tile boundary crosses a building and a road, so the discontinuity lands on high-contrast detail. In the second, a minimum-cost path follows low-contrast ground between them. STRAIGHT CUT The cut crosses a roof and a road; both show the radiometric step. ROUTED CUT A minimum-cost path avoids structures and crosses uniform ground. Routing the cut is free. Every other seam treatment costs either detail or contrast.

Core Algorithmic Pipeline

  1. Compute the overlap region between the incoming capture and the existing mosaic, in pixel space on the mosaic’s grid.
  2. Route a cut line through the overlap that prefers low-contrast ground — a minimum-cost path where cost is local image gradient.
  3. Histogram-match the incoming capture to the mosaic, using only pixels inside the overlap so the statistics describe the same ground.
  4. Feather across the cut line over a band defined in metres and converted to pixels.
  5. Verify along the seam by measuring the cross-line gradient against the local background gradient.
Match first, then feather — the order matters A chain of four stages: compute the overlap, route the cut line through low-contrast ground, histogram-match the incoming capture on overlap pixels only, and feather across a narrow band sized in metres. Overlap pixels valid in both region Cut line minimum-cost path mask Match fit on overlap only corrected Feather narrow band, in metres Feathering an unmatched radiometric gap produces a smudge instead of a seam, which is not an improvement.

Working Implementation

"""Seam-aware merge of an incoming orthophoto capture into an existing mosaic."""
from __future__ import annotations

import numpy as np
import rasterio
from rasterio.windows import from_bounds
from skimage.graph import route_through_array


def overlap_window(mosaic, incoming):
    """Pixel window, on the mosaic grid, where the two datasets both have data."""
    left = max(mosaic.bounds.left, incoming.bounds.left)
    bottom = max(mosaic.bounds.bottom, incoming.bounds.bottom)
    right = min(mosaic.bounds.right, incoming.bounds.right)
    top = min(mosaic.bounds.top, incoming.bounds.top)
    if right <= left or top <= bottom:
        raise ValueError("captures do not overlap; nothing to blend")
    return from_bounds(left, bottom, right, top, mosaic.transform)


def cut_line(a: np.ndarray, b: np.ndarray) -> np.ndarray:
    """Minimum-cost path down the overlap, preferring low-contrast ground.

    Cost combines how much the two captures disagree with how much local
    detail exists — a cut through a uniform field where both agree is free,
    a cut across a roof ridge is expensive.
    """
    disagreement = np.abs(a.astype("float32") - b.astype("float32")).mean(axis=0)
    gy, gx = np.gradient(a.mean(axis=0).astype("float32"))
    detail = np.hypot(gy, gx)
    cost = 1.0 + disagreement + 0.5 * detail

    rows, cols = cost.shape
    path, _ = route_through_array(
        cost, start=(0, cols // 2), end=(rows - 1, cols // 2), fully_connected=True
    )
    mask = np.zeros((rows, cols), dtype="uint8")
    for r, c in path:
        mask[r, :c] = 1                       # left of the cut belongs to `a`
    return mask


def match_histogram(source: np.ndarray, reference: np.ndarray,
                    valid: np.ndarray) -> np.ndarray:
    """Linear per-band match of source to reference, fitted on overlap pixels only."""
    matched = np.empty_like(source)
    for band in range(source.shape[0]):
        s = source[band][valid].astype("float64")
        r = reference[band][valid].astype("float64")
        if s.std() < 1e-6:
            matched[band] = source[band]
            continue
        gain = r.std() / s.std()
        offset = r.mean() - gain * s.mean()
        adjusted = source[band].astype("float64") * gain + offset
        info = np.iinfo(source.dtype)
        matched[band] = np.clip(adjusted, info.min, info.max).astype(source.dtype)
    return matched


def feather_mask(mask: np.ndarray, width_px: int) -> np.ndarray:
    """Turn a hard 0/1 mask into a smooth ramp of the requested width."""
    from scipy.ndimage import distance_transform_edt

    inside = distance_transform_edt(mask)
    outside = distance_transform_edt(1 - mask)
    signed = inside - outside
    return np.clip(0.5 + signed / (2.0 * max(width_px, 1)), 0.0, 1.0).astype("float32")


def merge_capture(mosaic_path, incoming_path, out_path, feather_m=6.0,
                  match=True):
    with rasterio.open(mosaic_path) as mosaic, rasterio.open(incoming_path) as incoming:
        window = overlap_window(mosaic, incoming)
        base = mosaic.read(window=window)
        new = incoming.read(window=from_bounds(*rasterio.windows.bounds(
            window, mosaic.transform), incoming.transform))

        valid = (base.max(axis=0) > 0) & (new.max(axis=0) > 0)
        if valid.sum() < 1000:
            raise ValueError("overlap holds too few valid pixels to fit a match")

        if match:
            new = match_histogram(new, base, valid)

        hard = cut_line(base, new)
        px_per_m = 1.0 / abs(mosaic.transform.a)
        alpha = feather_mask(hard, int(round(feather_m * px_per_m)))

        blended = (base.astype("float32") * alpha
                   + new.astype("float32") * (1.0 - alpha))

        profile = mosaic.profile
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(mosaic.read())
            dst.write(blended.astype(base.dtype), window=window)
    return out_path

Expressing the feather in metres and converting with the transform is what lets the same call work on a 5 cm drone capture and a 50 cm aerial survey. A width hard-coded in pixels produces a two-metre smear on one and an invisible ramp on the other.

Validation & Output Verification

A seam is a local defect, so the check has to be local. Whole-image statistics will not see a one-pixel step:

import numpy as np
import rasterio

def seam_gradient_ratio(path, mask, band=1, halo=24):
    """Mean cross-seam gradient divided by mean local gradient.

    Close to 1.0 means the seam is indistinguishable from ordinary image
    texture. Above ~1.5 it is visible on screen.
    """
    with rasterio.open(path) as src:
        arr = src.read(band).astype("float32")

    edge = np.abs(np.gradient(mask.astype("float32"))[1]) > 0
    gy, gx = np.gradient(arr)
    grad = np.hypot(gy, gx)

    from scipy.ndimage import binary_dilation
    near = binary_dilation(edge, iterations=2)
    background = binary_dilation(edge, iterations=halo) & ~near

    return float(grad[near].mean() / max(grad[background].mean(), 1e-6))
# The merged mosaic must still be a valid COG with its pyramid intact
python -m rio_cogeo validate merged_mosaic.tif

# And the seam must not be the brightest edge in its neighbourhood
python -c "
from seam import seam_gradient_ratio
import numpy as np
mask = np.load('cut_mask.npy')
r = seam_gradient_ratio('merged_mosaic.tif', mask)
print(f'seam gradient ratio {r:.2f}')
assert r < 1.5, 'seam is visible — widen the feather or re-run histogram matching'
"

Record the ratio in the merge manifest alongside the arbitration decisions described in the parent guide. It turns “the mosaic looks fine” into a number that a later merge can be compared against.

Choosing the feather width Horizontal bars showing the measured seam gradient ratio at four feather widths: a hard cut is clearly visible, six metres brings the ratio close to background texture, and wide feathers reduce it no further while smearing real edges. SEAM GRADIENT RATIO BY FEATHER WIDTH 0 m (hard cut) 3.4 clearly visible on screen 2 m 1.9 6 m 1.2 indistinguishable from image texture 20 m 1.1 no further gain, and road markings blur The curve flattens quickly, so anything past the third bar costs detail and buys nothing.

Failure Modes

  • The blended band is a visible smudge — symptom: no hard line, but an obvious soft stripe. Root cause: feathering a large radiometric difference without matching first. Fix: histogram-match on overlap pixels, then feather narrowly.

  • Real detail lost along the seam — symptom: road markings fade near the boundary. Root cause: feather width far exceeds feature size. Fix: express the width in metres, keep it under about ten metres, and route the cut line away from detailed ground.

  • The cut line runs straight through a building — symptom: half a roof from each capture, with an offset ridge. Root cause: cost function weighted disagreement only. Fix: include the local gradient term so structures are expensive to cross.

  • Histogram match makes the incoming tile obviously wrong — symptom: the whole tile shifts colour. Root cause: statistics fitted over the full tile rather than the overlap, so different ground cover drove the fit. Fix: restrict the fit to pixels valid in both captures, as the implementation does.

Back to Merge Strategies for Raster and Imagery Mosaics