Correcting Vignetting and Lens Falloff

The reflectance orthomosaic has a regular pattern of darker patches across it. They are not clouds — their edges are straight, they repeat at a fixed spacing, and the spacing matches the distance between frame centres rather than anything on the ground. In the index raster derived from it, the same pattern appears as a systematic bias that shifts a plot’s mean depending on where in a frame that plot happened to fall.

This is lens falloff, and on the small optics multispectral cameras carry it reaches 25–40 % at the frame corners. Uncorrected it is one of the largest errors in the whole radiometric chain — larger than most of the effects the chain is carefully removing.

# Where the falloff comes from, and why it is not uniform

Two mechanisms combine. Optical vignetting is physical obstruction: light entering at an angle is partly blocked by the lens barrel, and the effect grows with the angle, hence with distance from the optical axis. Natural falloff is geometric — the well-known cosine-fourth law — and applies to any lens regardless of construction.

Both are radial about the optical centre, which is not the centre of the sensor. Manufacturing tolerances put it tens of pixels away, and correcting about the image centre instead leaves a residual gradient across the frame. That residual is small — a few percent — but it is systematic and it does not average out across a mosaic.

A third contribution is band-dependent. Each band on a multispectral rig has its own lens, so each has its own falloff profile, and applying one band’s correction to another introduces a spectral error that varies with position in the frame. That is the worst kind: it looks like a real spatial pattern in the index.

Figure 1 — Five bands, five lenses, five profiles. One correction for all of them is not an approximation; it is a new error.

# Minimal reproducible solution

Where the vendor supplies coefficients, use them — they were measured on that unit. Where they do not, fit from a flat field.

import numpy as np


def vignette_map(shape: tuple[int, int], *, centre: tuple[float, float],
                 coeffs: list[float]) -> np.ndarray:
    """Per-pixel multiplicative correction from a radial polynomial.

    The polynomial expresses the response as a function of radius, so the
    correction is its reciprocal. Clipping the denominator keeps a badly
    conditioned polynomial from producing infinities in the far corners.
    """
    h, w = shape
    cy, cx = centre
    yy, xx = np.mgrid[0:h, 0:w]
    r = np.hypot(xx - cx, yy - cy)

    response = np.zeros_like(r, dtype=np.float64)
    for power, c in enumerate(coeffs):
        response += c * r ** power

    response = np.clip(response, 0.2, 1.5)
    return 1.0 / response


def fit_vignette(flat: np.ndarray, *, order: int = 4,
                 centre: tuple[float, float] | None = None) -> dict:
    """Fit a radial response model to a flat-field image.

    Fitting a smooth model rather than dividing by the flat field directly is
    what keeps the flat field's own noise and dust specks out of every
    calibrated frame in the survey.
    """
    h, w = flat.shape
    cy, cx = centre if centre else (h / 2.0, w / 2.0)
    yy, xx = np.mgrid[0:h, 0:w]
    r = np.hypot(xx - cx, yy - cy).ravel()
    v = flat.ravel().astype(np.float64)

    good = np.isfinite(v) & (v > 0)
    design = np.vander(r[good], order + 1, increasing=True)
    coeffs, residuals, *_ = np.linalg.lstsq(design, v[good], rcond=None)
    coeffs = coeffs / coeffs[0]              # normalise so the centre reads 1.0

    model = (np.vander(r, order + 1, increasing=True) @ coeffs).reshape(h, w)
    rms = float(np.sqrt(np.mean((flat[good.reshape(h, w)] -
                                 model[good.reshape(h, w)] * flat.max()) ** 2)))
    return {"coeffs": coeffs.tolist(), "centre": (cy, cx), "fit_rms": rms}

Normalising the coefficients so the centre reads exactly 1.0 matters for the same reason the panel step does: the vignette correction must change the shape of the frame’s response without changing its overall level, or it silently rescales every reflectance value alongside the panel factor.

# Finding the optical centre

Correcting about the wrong centre leaves a gradient. The centre can be recovered from the flat field by minimising the residual of the radial fit over candidate centres.

import numpy as np
from scipy.optimize import minimize


def find_optical_centre(flat: np.ndarray, order: int = 4) -> tuple[float, float]:
    """Locate the centre that makes the response most nearly radial.

    A response that is genuinely radial about some point will fit a radial
    polynomial with a small residual only when that point is used. Searching
    over the centre is therefore a well-posed problem with one clear minimum.
    """
    h, w = flat.shape

    def cost(params):
        cy, cx = params
        yy, xx = np.mgrid[0:h, 0:w]
        r = np.hypot(xx - cx, yy - cy).ravel()
        v = flat.ravel().astype(np.float64)
        good = np.isfinite(v) & (v > 0)
        design = np.vander(r[good], order + 1, increasing=True)
        coeffs, *_ = np.linalg.lstsq(design, v[good], rcond=None)
        return float(np.mean((design @ coeffs - v[good]) ** 2))

    res = minimize(cost, x0=np.array([h / 2, w / 2]), method="Nelder-Mead",
                   options={"xatol": 0.5, "fatol": 1e-3})
    return float(res.x[0]), float(res.x[1])

On a typical multispectral sensor this lands 10–60 pixels from the image centre, and correcting about it removes a residual gradient of one to three percent — small, but exactly the size of the effects the survey is measuring.

Figure 3 — What a correct falloff profile does to the radial trace.

# Edge-case matrix

# Verification snippet

The definitive check is that the same ground looks the same regardless of where in a frame it fell.

import numpy as np


def position_bias(values: np.ndarray, frame_x: np.ndarray, frame_y: np.ndarray,
                  centre: tuple[float, float], bins: int = 5) -> dict:
    """Does a calibrated value depend on where in the frame it was observed?

    Bins observations of the same ground by their radial position within the
    source frame. After a correct vignette correction the bins agree; before
    it, the outer bins read systematically low.
    """
    cy, cx = centre
    r = np.hypot(frame_x - cx, frame_y - cy)
    edges = np.percentile(r, np.linspace(0, 100, bins + 1))

    means = []
    for i in range(bins):
        sel = (r >= edges[i]) & (r <= edges[i + 1]) & np.isfinite(values)
        means.append(float(np.median(values[sel])) if sel.sum() > 50 else np.nan)

    means_arr = np.asarray(means)
    trend = float(np.nanmax(means_arr) - np.nanmin(means_arr))
    baseline = float(np.nanmedian(means_arr))
    return {"bin_medians": means, "relative_trend": trend / max(baseline, 1e-9),
            "ok": trend / max(baseline, 1e-9) < 0.02}

A relative trend under about 2 % across the radial bins means the correction is working. Anything larger and the residual is still in the data, and it will appear in the mosaic wherever frame edges overlap frame centres.

Figure 2 — The definitive test, and the one that needs no external reference.

# Producing a usable flat field in the field

Where vendor coefficients are unavailable, the flat field has to be captured, and it is easier than it sounds. The requirement is a uniformly lit, uniformly reflective surface filling the frame — which a large reflectance panel, an overcast sky, or a sheet of white card under diffuse light all provide.

Two practical points decide whether the result is usable. Defocus slightly, so that any texture in the target is blurred below the scale of the model being fitted; the model is smooth and radial, so fine detail in the target only adds noise to the fit. And capture several frames and average them, which removes sensor noise without any additional processing.

An overcast sky is the most accessible option and works well: point the camera up on a uniformly grey day, well away from the sun’s position, and the illumination across the frame is flat to within a percent or so. Avoid clear blue sky, which has a strong gradient with angle from the sun and will be fitted as though it were falloff.

# When to escalate

• The fit residual is large and the target was genuinely flat. The response may not be radial — a decentred or damaged lens produces an asymmetric pattern that no radial model captures. A two-dimensional surface fit will describe it, but the underlying fault is hardware.
• Correction improves the mosaic and the index still shows frame patterns. Check band alignment next; misregistration produces its own frame-correlated artefacts, covered in fixing band misregistration artifacts in index rasters.
• The pattern changed between two flights with the same camera. The aperture or a lens element moved. Re-derive the coefficients rather than reusing them.

# Related

• Radiometric Calibration in Python
• Fixing Exposure Drift Across a Multispectral Flight
• Fixing Visible Seamlines and Colour Banding

← Radiometric Calibration in Python