Co-registering Bands with OpenCV Homographies

The obvious approach to aligning two bands — cross-correlate them and find the offset — works poorly on multispectral imagery, and the reason is instructive. Vegetation is dark in red and bright in near-infrared. Soil is the other way round. A correlation between the two bands over a field of crop and bare patches is negative, and a matcher looking for similarity finds the alignment that maximises anti-correlation, which is not the alignment that puts the same ground in the same place.

What survives a wavelength change is structure: the position of an edge is identical in every band even when the sign of the contrast is not. This page builds an alignment on that property using OpenCV, and — more importantly — measures whether the result is good enough to use. It is the implementation detail behind band alignment and stacking for multispectral sets.

Why gradient features work where intensity does not

A feature detector such as AKAZE or ORB responds to local structure — corners, blobs, edges — detected in a scale space built from gradients. The descriptor it attaches encodes the pattern of gradients around the keypoint, not their absolute values. When the contrast of a region reverses between bands, the gradient magnitudes are preserved and the signs flip, and a binary descriptor built from comparisons is affected far less than an intensity template would be.

That is why these detectors work across wavelengths at all, and it also explains their failure mode: they fail where there is no structure. A frame over uniform water, a bare tilled field at low sun, or a dense uniform canopy presents few reliable keypoints in any band, and the match count collapses.

The practical response is to normalise each band independently before detection — stretching each to the same nominal range so the detector’s internal thresholds behave the same way on both — and to treat a low match count as a signal to fall back rather than as something to force.

Intensity reverses between bands; edge position does not Two profiles across the same boundary between crop and soil. In the red band the crop is dark and the soil bright, so the profile steps up at the boundary. In the near-infrared band the crop is bright and the soil dark, so the profile steps down at the same position. Beneath, the gradient magnitude of both profiles peaks at exactly the same place. A note states that a correlation between the two bands is negative while their gradient magnitudes agree, which is why gradient-based features work across wavelengths. boundary red — crop dark, soil bright near-infrared — reversed gradient magnitude — identical in both Correlation between the bands is negative here. Their structure agrees exactly.

Figure 1 — The property that makes cross-band matching possible at all.

Minimal reproducible solution

import cv2
import numpy as np


def coregister(reference: np.ndarray, target: np.ndarray, *,
               ransac_px: float = 2.0, ratio: float = 0.75) -> dict:
    """Estimate the homography mapping `target` onto `reference`.

    Each band is stretched to the same nominal range before detection so the
    detector's internal contrast thresholds behave identically on both — a
    step that roughly doubles the usable match count on a typical frame pair.
    """
    def prep(img: np.ndarray) -> np.ndarray:
        a = img.astype(np.float32)
        lo, hi = np.percentile(a[np.isfinite(a)], [2, 98])
        return np.clip((a - lo) / max(hi - lo, 1e-6) * 255, 0, 255).astype(np.uint8)

    detector = cv2.AKAZE_create(threshold=0.0008)
    ka, da = detector.detectAndCompute(prep(reference), None)
    kb, db = detector.detectAndCompute(prep(target), None)
    if da is None or db is None:
        raise ValueError("no descriptors in one of the bands")

    matcher = cv2.BFMatcher(cv2.NORM_HAMMING)
    pairs = matcher.knnMatch(db, da, k=2)
    good = [m for m, n in pairs if m.distance < ratio * n.distance]
    if len(good) < 25:
        raise ValueError(f"only {len(good)} matches survive the ratio test")

    src = np.float32([kb[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)
    dst = np.float32([ka[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)

    H, mask = cv2.findHomography(src, dst, cv2.RANSAC, ransac_px,
                                 maxIters=5000, confidence=0.999)
    if H is None:
        raise ValueError("RANSAC found no consistent homography")

    keep = mask.ravel().astype(bool)
    resid = np.linalg.norm(
        cv2.perspectiveTransform(src[keep], H) - dst[keep], axis=2).ravel()

    return {"H": H, "matches": len(good), "inliers": int(keep.sum()),
            "inlier_fraction": float(keep.mean()),
            "median_residual_px": float(np.median(resid)),
            "p95_residual_px": float(np.percentile(resid, 95))}

Choosing the RANSAC threshold

The threshold decides which correspondences count as agreeing with the model, and getting it wrong fails in both directions. Too tight and RANSAC discards good matches, leaving the fit under-constrained and unstable between frames. Too loose and it admits mismatches, which drag the homography and produce a transform that is smooth, plausible and wrong.

A useful discipline is to sweep it once on a representative frame pair and look for the plateau: over a range of thresholds the inlier count rises steadily and the residual stays flat, and the right value sits at the top of that plateau.

import numpy as np


def sweep_ransac_threshold(src: np.ndarray, dst: np.ndarray,
                           thresholds=(0.75, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0)) -> list[dict]:
    """Inliers and residual against the RANSAC threshold, to find the plateau."""
    import cv2
    rows = []
    for t in thresholds:
        H, mask = cv2.findHomography(src, dst, cv2.RANSAC, t, maxIters=5000)
        if H is None:
            rows.append({"threshold": t, "inliers": 0, "median_residual_px": None})
            continue
        keep = mask.ravel().astype(bool)
        resid = np.linalg.norm(
            cv2.perspectiveTransform(src[keep], H) - dst[keep], axis=2).ravel()
        rows.append({"threshold": t, "inliers": int(keep.sum()),
                     "inlier_fraction": float(keep.mean()),
                     "median_residual_px": float(np.median(resid))})
    return rows
Inlier count and residual against the RANSAC threshold Two curves against a RANSAC threshold from zero point seven five to six pixels. The inlier count rises steeply from about forty at the tightest threshold to a plateau of around two hundred and forty between one and a half and three pixels, then continues rising as mismatches are admitted. The median residual stays flat at about a third of a pixel across the plateau and then rises sharply beyond three pixels. A shaded band marks the plateau as the region where the threshold should be chosen. plateau — choose here 0.75 1.0 1.5 2.0 3.0 6.0 px RANSAC inlier threshold inliers median residual Past the plateau the inlier count keeps rising because mismatches are being admitted.

Figure 2 — Reading the plateau. A rising inlier count with a rising residual is not an improvement.

A homography that is correct and one that has merely converged Two columns. The correct column notes that the transform is close to a small translation with a slight scale, that residuals are distributed evenly across the frame, that it is consistent frame to frame across the flight, and that applying it visibly removes boundary artefacts from the index. The merely converged column notes a large rotation or shear that no rigid band offset could produce, residuals concentrated in one region, a transform that differs wildly between consecutive frames, and an index whose boundary artefacts move rather than disappear. correct a small translation, slight scale residuals even across the frame consistent frame to frame boundary artefacts disappear merely converged large rotation or shear residuals concentrated in one region wildly different between frames artefacts move rather than vanish The band offset is a physical property of the sensor: it barely changes between frames.

Figure 3 — Frame-to-frame consistency is the cheapest validity check available.

Edge-case matrix

Situation Symptom Handling
Uniform water or canopy Few keypoints, matching fails Phase correlation fallback
Motion blur during a turn Descriptors unstable Reject the frame
Strong relief in frame Good global fit, bad over structures Homography cannot model parallax
Ratio test too strict Under 25 matches Relax to 0.8, or change detector
Ratio test too loose Many matches, low inlier fraction Tighten; mismatches drag the fit
Bands at very different resolution Scale mismatch Resample to a common grid first
Repeated structures (crop rows) Confident wrong matches Constrain with the expected offset range
Frame near the survey edge Half the frame has no overlap Expected; judge on inlier count, not fraction

The repeated-structures row is worth attention on agricultural sites. Crop rows are a regular pattern, and a matcher can lock onto the wrong row with high confidence. Constraining the search to transforms near the expected sensor offset — which is known from the rig geometry to within a pixel or two — eliminates the failure entirely.

Verification snippet

import numpy as np


def sanity_check_homography(H: np.ndarray, shape: tuple[int, int],
                            *, max_shift_px: float = 60.0,
                            max_scale_deviation: float = 0.05) -> list[str]:
    """Is this transform physically plausible for two bands of one rig?

    Bands on a fixed rig differ by a small translation and a tiny scale
    change. A homography implying a large rotation or a ten percent scale is
    a fitting artefact, however good its residual looks.
    """
    h, w = shape
    problems = []

    corners = np.float32([[0, 0], [w, 0], [w, h], [0, h]]).reshape(-1, 1, 2)
    import cv2
    moved = cv2.perspectiveTransform(corners, H).reshape(-1, 2)
    shifts = np.linalg.norm(moved - corners.reshape(-1, 2), axis=1)
    if shifts.max() > max_shift_px:
        problems.append(f"corner moves {shifts.max():.0f} px — implausible for a fixed rig")

    scale = np.sqrt(abs(np.linalg.det(H[:2, :2])))
    if abs(scale - 1.0) > max_scale_deviation:
        problems.append(f"implied scale {scale:.3f} — bands should be near unity")

    shear = abs(H[0, 1]) + abs(H[1, 0])
    if shear > 0.05:
        problems.append(f"implied rotation or shear of {shear:.3f}")
    return problems

Checking physical plausibility alongside the residual catches the case where RANSAC has found an internally consistent set of mismatches. A transform that fits its inliers beautifully and moves the frame corner by two hundred pixels is wrong regardless of how good the residual is.

Running it over a whole flight efficiently

Estimating a homography per band pair per frame is four estimations for a five-band rig, times a few thousand frames. Two adjustments keep that to a few minutes rather than an hour.

Detect keypoints in the reference band once per frame and reuse them for all four pairs, rather than re-detecting for each. Detection dominates the runtime, and the reference band’s keypoints do not change between pairs.

Seed each frame’s estimate from the previous frame’s transform and reject solutions far from it. Consecutive frames on a line differ by a fraction of a pixel, so the seed is almost the answer, and the rejection doubles as the repeated-structure guard described above.

When to escalate

  • Matching fails on a large share of frames. The imagery is the problem — motion blur from flying too fast for the exposure, or a site with genuinely no structure. Neither is fixed in processing.
  • Residuals are acceptable and index fringes persist. Look at whether alignment was applied before or after orthorectification, and whether the bands were resampled a different number of times.
  • The rig’s geometry appears to have changed. A dropped or serviced camera can shift a sensor. Compare the median transform across a flight with a previous survey’s; a consistent change is a hardware event worth recording.

← Band Alignment and Stacking for Multispectral Sets