Fixing Band Misregistration Artifacts in Index Rasters

The NDVI raster has a bright line along the north side of every hedge and a dark line along the south side. Around each bare patch in the field there is a thin rim of implausibly high values. The agronomist asks what the bright lines mean, and the honest answer is that they do not mean anything: they are an arithmetic consequence of two bands that do not quite overlap.

The artefact is easy to recognise once seen, easy to measure, and impossible to remove from the index raster itself. The fix is always upstream, in the alignment covered by band alignment and stacking for multispectral sets.

# The arithmetic that produces the fringe

Consider a boundary between crop and soil, and suppose the near-infrared band is shifted two pixels north relative to red. At the boundary there is a strip two pixels wide where the near-infrared value still belongs to the crop while the red value already belongs to the soil.

For a normalised difference index, that strip combines the highest near-infrared with the highest red-contrast — crop NIR over soil red — and the index there exceeds anything real in the scene. Two pixels away on the other side, the combination is reversed and the index drops below anything real.

Two properties follow, and together they identify the artefact conclusively. The fringe is always paired: a bright line on one side of an edge and a dark line on the other, with the pairing oriented along the misregistration direction. And it is oriented consistently across the whole raster, because the offset between two sensors is the same everywhere in a frame.

Figure 1 — The arithmetic, and the signature it leaves in plan view.

# Minimal reproducible solution

The offset can be measured from the index raster itself, which is useful when the source frames are no longer to hand: the direction in which index extremes concentrate at edges is the misregistration direction.

import numpy as np
from scipy import ndimage


def estimate_misregistration(index: np.ndarray, *, edge_percentile: float = 92.0,
                             search_px: int = 4) -> dict:
    """Recover the band offset direction and magnitude from the index raster.

    Extreme index values cluster on one side of every edge. Correlating the
    extremes against a directional gradient recovers the offset without the
    source bands.
    """
    a = np.nan_to_num(index.astype(np.float32))
    gy, gx = np.gradient(a)
    grad = np.hypot(gx, gy)
    edges = grad > np.nanpercentile(grad, edge_percentile)
    if edges.sum() < 500:
        return {"note": "not enough edge pixels to judge"}

    extreme_hi = a > np.nanpercentile(a[np.isfinite(a)], 99.5)
    best = None
    for dy in range(-search_px, search_px + 1):
        for dx in range(-search_px, search_px + 1):
            if dy == 0 and dx == 0:
                continue
            shifted = ndimage.shift(edges.astype(np.float32), (dy, dx), order=0)
            overlap = float(np.mean(shifted[extreme_hi]))
            if best is None or overlap > best[0]:
                best = (overlap, dy, dx)

    overlap, dy, dx = best
    return {"offset_px": (dy, dx),
            "magnitude_px": float(np.hypot(dy, dx)),
            "confidence": overlap,
            "note": "offset direction recovered from the index raster alone"}

The result identifies the problem and does not fix it. Which brings us to the part that matters.

# Why the index raster cannot be repaired

It is tempting to filter the fringes out — a median filter, or masking the extreme percentiles, or eroding the plot boundaries away from edges. None of these is a fix, for a reason worth being precise about.

The fringe is not noise added to a correct value. It is what you get when two different places on the ground are divided by each other. The information needed to produce the correct index at those pixels — the near-infrared value at the same ground position as the red value — is not present in the raster at all; it is in the source frames, two pixels away.

Filtering therefore replaces a wrong value with a smoothed wrong value. It looks better, the extremes disappear, and the plot statistics are still biased because the bias was never in the extremes alone. The full-width effect of a two-pixel offset extends across every transition zone in the scene, not merely the pixels that clipped.

import numpy as np


def fringe_contribution(index: np.ndarray, edges: np.ndarray,
                        offset_px: float) -> dict:
    """How much of a plot's mean comes from the affected transition zone.

    The affected width is the misregistration magnitude, so on a plot with
    many internal edges — crop rows, tramlines — the affected fraction can
    be large even for a sub-pixel offset.
    """
    from scipy import ndimage
    affected = ndimage.binary_dilation(edges, iterations=max(int(round(offset_px)), 1))
    valid = np.isfinite(index)
    frac = float(np.count_nonzero(affected & valid) / max(np.count_nonzero(valid), 1))

    clean_mean = float(np.nanmean(index[valid & ~affected]))
    all_mean = float(np.nanmean(index[valid]))
    return {"affected_fraction": frac,
            "mean_all": all_mean, "mean_clean": clean_mean,
            "bias": all_mean - clean_mean}

Running that on a real plot is often the argument that settles the discussion: on a field with tramlines every twenty-four metres, a one-pixel offset at 5 cm resolution can affect several percent of the area and shift the plot mean by more than the seasonal change being measured.

Figure 3 — Four questions, and the last two matter most on undulating ground.

# Edge-case matrix

The last row is the trap. A pipeline that masks extreme index values before anyone looks at the raster hides the diagnostic while keeping the error, which is the worst of both.

# Verification snippet

import numpy as np
from scipy import ndimage


def fringe_test(index: np.ndarray, *, edge_percentile: float = 92.0) -> dict:
    """Do extreme index values cluster on one side of edges?

    A real feature is not systematically one-sided. A misregistration fringe
    is, and the asymmetry between the two sides of an edge is a direct,
    single-number test.
    """
    a = np.nan_to_num(index.astype(np.float32))
    gy, gx = np.gradient(a)
    grad = np.hypot(gx, gy)
    edges = grad > np.nanpercentile(grad, edge_percentile)

    hi = a > np.nanpercentile(a[np.isfinite(a)], 99.0)
    lo = a < np.nanpercentile(a[np.isfinite(a)], 1.0)

    near = ndimage.binary_dilation(edges, iterations=2)
    hi_near = float(np.count_nonzero(hi & near) / max(np.count_nonzero(hi), 1))
    lo_near = float(np.count_nonzero(lo & near) / max(np.count_nonzero(lo), 1))

    return {"extreme_high_near_edges": hi_near,
            "extreme_low_near_edges": lo_near,
            "misregistration_likely": hi_near > 0.6 and lo_near > 0.6,
            "note": ("extremes concentrate at edges — band misregistration"
                     if hi_near > 0.6 and lo_near > 0.6
                     else "extremes are distributed; edges are probably real")}

Running this as a gate on every index raster, before any masking, catches the problem at the point where re-aligning is still a rerun rather than a re-flight.

Figure 2 — Why “it is only a pixel” depends entirely on how many edges the field has.

# When to escalate

• The source frames are gone. The index cannot be corrected. Report the measured offset and its estimated bias alongside the product, and re-process from frames next time.
• Alignment is good and fringes persist. Check whether the index was computed from an orthomosaic assembled from differently aligned frames — a per-frame alignment applied after mosaicking does not help.
• The offset changed partway through the flight. A sensor has moved physically. That is a maintenance finding, and every subsequent survey with that rig is affected until it is fixed.

# Related

• Band Alignment and Stacking for Multispectral Sets
• Co-registering Bands with OpenCV Homographies
• Computing Vegetation Indices with Rasterio

← Band Alignment and Stacking for Multispectral Sets