noaa-goes-visualization/suvi/detectors.py

"""Bad-frame detectors.

Every detector is a pure function: data in, a :class:`Verdict` out.  Nothing here
renames a file, writes a plot, or touches the database -- that is the driver's job.
This is what lets the same code run in production (``filter_FITS.py``) and under the
test bench, and lets several detectors judge the same frame independently.

Two shapes of detector:

* **Frame-level** (:func:`header_v1`, :func:`geometry_v1`) judge one frame alone.
* **Series-level** (:func:`temporal_v1`, :func:`crosssat_v1`) judge a frame in the
  context of its neighbours in time, or of the other satellite at the same instant.
  They work from :class:`FrameFeatures` -- cached header values plus a small
  thumbnail -- rather than full 1280x1280 arrays, which keeps a whole day of context
  in memory at once.

Every detector reports continuous ``scores`` alongside its verdict.  The driver
stores those, so re-tuning a threshold later is a database query rather than another
pass over the archive.
"""

from dataclasses import dataclass, field

import cv2 as cv
import numpy as np

# --------------------------------------------------------------------------- types


@dataclass(frozen=True)
class Verdict:
    """One detector's judgement of one frame."""

    verdict: str  # 'good' | 'bad' | 'unknown'
    reason: str | None = None
    scores: dict = field(default_factory=dict)


GOOD = "good"
BAD = "bad"
UNKNOWN = "unknown"


@dataclass
class FrameFeatures:
    """The cheap summary of a frame that series-level detectors work from."""

    slot: tuple  # (satellite, wavelength, t_start)
    header: dict = field(default_factory=dict)
    #: Read error from fitsio, if the frame could not be parsed at all.
    error: str | None = None
    #: Small float32 image (see THUMBNAIL_SIZE), or None if not loaded.
    thumbnail: np.ndarray | None = None

    @property
    def t_start(self):
        return self.slot[2]

    @property
    def wavelength(self):
        return self.slot[1]


#: Series detectors compare thumbnails, not full frames.  128px preserves the disc,
#: active regions and gross structure while making a day of context cheap to hold.
THUMBNAIL_SIZE = 128


def thumbnail(image, size=THUMBNAIL_SIZE):
    """Downsample a frame for structural comparison.  NaNs become zero."""
    clean = np.nan_to_num(np.asarray(image, dtype=np.float32), nan=0.0, posinf=0.0, neginf=0.0)
    return cv.resize(clean, (size, size), interpolation=cv.INTER_AREA)


# ------------------------------------------------------------------ header_v1

#: Plausible radiance range of a healthy frame, per passband, in W m-2 sr-1.
#:
#: Deliberately wide -- roughly a factor of five outside the observed annual range of
#: good frames.  Sampled across 24 days of 2024 on both satellites, the *rejected*
#: frames' IMG_MEAN distribution sits almost entirely inside the good one (171A: good
#: 0.288-0.452, rejected median 0.352), so a tight bound here would cost enormous
#: numbers of false positives while catching little.  What this rule is actually for
#: is gross dropout and blowout -- eclipse frames run 1e-4, three orders of magnitude
#: below the floor -- and it is priced accordingly.
HEADER_MEAN_BOUNDS = {
    94: (0.004, 0.41),
    131: (0.005, 0.39),
    171: (0.057, 2.3),
    195: (0.15, 5.4),
    284: (0.13, 7.0),
    304: (0.29, 20.0),
}

#: Solar diameter in pixels; varies ~3% over the year with Earth-Sun distance.
DIAM_SUN_RANGE = (740.0, 800.0)
#: How far the recorded sun centre may sit from the image centre, in pixels.
MAX_CRPIX_OFFSET = 8.0


def header_v1(features, mean_bounds=None, diam_range=DIAM_SUN_RANGE,
              max_crpix_offset=MAX_CRPIX_OFFSET):
    """Judge a frame from its header alone -- no pixel decompression.

    The instrument already reports most of what matters: ``ECLIPSE`` marks frames
    taken through Earth's shadow, ``EMPTY`` marks frames built from no source
    imagery, and ``IMG_MEAN``/``IMG_SDEV`` summarise the radiance.  On sampled
    archive days every frame the pixel filter rejected also carried ``ECLIPSE=2``,
    at roughly a fifth of the cost of decompressing the image.

    ``DEGRADED`` is recorded as a score but is deliberately *not* a rejection rule.
    It is set for reasons unrelated to a frame's usability: across sampled days the
    195A band carried ``DEGRADED=T`` on 40/40 good frames on 2024-01-20 and 0/40 on
    2024-07-04, on both satellites.  Rejecting on it would discard entire months of
    that band.  The bench can re-weight it from the stored score if that turns out
    to be worth doing.
    """
    mean_bounds = mean_bounds or HEADER_MEAN_BOUNDS
    header = features.header
    scores = {}

    if features.error:
        return Verdict(BAD, f"unreadable: {features.error}", {"read_ok": 0.0})
    if not header:
        return Verdict(UNKNOWN, "no header cached", scores)

    for field_name in ("img_mean", "img_sdev", "diam_sun", "eclipse", "degraded"):
        if field_name in header:
            scores[field_name] = float(header[field_name])

    if header.get("empty"):
        return Verdict(BAD, "EMPTY flag set", scores)
    eclipse = header.get("eclipse")
    if eclipse:
        return Verdict(BAD, f"ECLIPSE flag set ({int(eclipse)})", scores)

    # A frame built from no source images carries no information regardless of what
    # the radiance statistics happen to say.
    if header.get("num_imgs") is not None and header["num_imgs"] < 1:
        scores["num_imgs"] = float(header["num_imgs"])
        return Verdict(BAD, "no source images", scores)

    # The band the file claims must match the band its name puts it in, or the
    # archive has a misfiled frame and every downstream threshold is wrong.
    declared = header.get("wavelnth")
    if declared is not None and int(declared) != features.wavelength:
        scores["wavelnth"] = float(declared)
        return Verdict(BAD, f"wavelength {int(declared)} != {features.wavelength}", scores)

    mean = header.get("img_mean")
    if mean is not None:
        low, high = mean_bounds.get(features.wavelength, (0.0, np.inf))
        if mean < low:
            return Verdict(BAD, f"IMG_MEAN {mean:.3g} below {low:g}", scores)
        if mean > high:
            return Verdict(BAD, f"IMG_MEAN {mean:.3g} above {high:g}", scores)

    sdev = header.get("img_sdev")
    if sdev is not None and sdev <= 0.0:
        return Verdict(BAD, "zero variance", scores)

    diameter = header.get("diam_sun")
    if diameter is not None and not (diam_range[0] <= diameter <= diam_range[1]):
        return Verdict(BAD, f"DIAM_SUN {diameter:.1f} outside {diam_range}", scores)

    crpix1, crpix2 = header.get("crpix1"), header.get("crpix2")
    if crpix1 is not None and crpix2 is not None:
        offset = max(abs(crpix1 - 640.5), abs(crpix2 - 640.5))
        scores["crpix_offset"] = float(offset)
        if offset > max_crpix_offset:
            return Verdict(BAD, f"sun centre offset {offset:.1f}px", scores)

    return Verdict(GOOD, None, scores)


# ---------------------------------------------------------------- geometry_v1

#: Per-band radiance ceiling used to normalise the image before shape analysis.
#: Index order matches paths.WAVELENGTHS.
GEOMETRY_THRESHOLDS = {94: 0.050, 131: 0.10, 171: 1.00, 195: 1.40, 284: 1.0, 304: 2.50}

EXPECTED_DIMS = 1280
HALF_DIMS = EXPECTED_DIMS // 2
VALID_RADII = (383, 394)
AVG_RADIUS = (VALID_RADII[0] + VALID_RADII[1]) // 2
MAX_RADIUS_ERROR = 80
MAX_CENTER_SKEW = 7
RATIO_ABOVE_THRESH_MAX = 0.4
RATIO_ABOVE_THRESH_MIN = 0.07
EDGE_THRESH = 0.98
MAX_GOF = 0.2

_ideal_axis_cache = {}


def ideal_disc_axis():
    """Column-mean profile of a perfect solar disc, used as the shape reference."""
    if "axis" not in _ideal_axis_cache:
        disc = np.zeros((HALF_DIMS, HALF_DIMS))
        cv.circle(disc, (HALF_DIMS // 2, HALF_DIMS // 2), AVG_RADIUS // 2, 1, -1)
        _ideal_axis_cache["axis"] = np.average(disc, 0)
    return _ideal_axis_cache["axis"]


def geometry_v1(image, wavelength):
    """The pre-existing pixel-geometry filter, preserved as the baseline.

    Normalises against a per-band radiance ceiling, then checks the fraction of
    saturated pixels, the intensity-weighted centroid, the radius implied by the
    98th-percentile cumulative edges, and the fit against an ideal disc profile.

    Deliberately kept behaviourally identical to the original so bench numbers
    describe what the pipeline actually did, with two known weaknesses left intact
    for the bench to quantify:

    * the centre check is one-sided -- ``HALF_DIMS//2 - centre > skew`` catches a
      disc displaced up/left but not one displaced down/right;
    * non-finite pixels propagate as NaN through every comparison, and ``NaN > x``
      is False, so an all-NaN frame passes every check.

    The one deliberate change: the original raised on an unexpected image size and
    left the file unlabelled.  Here that is a ``bad`` verdict, because "no verdict"
    is not a usable baseline to score against.
    """
    scores = {}
    if image is None:
        return Verdict(BAD, "no image data", scores)
    if image.shape[0] != EXPECTED_DIMS or image.shape[1] != EXPECTED_DIMS:
        return Verdict(BAD, f"unexpected dimensions {image.shape}", scores)

    threshold = GEOMETRY_THRESHOLDS.get(wavelength)
    if threshold is None:
        return Verdict(UNKNOWN, f"no threshold for band {wavelength}", scores)

    data = cv.resize(
        np.asarray(image, dtype=np.float32),
        dsize=(HALF_DIMS, HALF_DIMS),
        interpolation=cv.INTER_LINEAR,
    )
    above = data > threshold
    normalised = np.copy(data)
    normalised[above] = threshold
    normalised /= threshold

    ratio = np.count_nonzero(above) / data.shape[0] / data.shape[1]
    scores["ratio_above_thresh"] = float(ratio)
    if ratio > RATIO_ABOVE_THRESH_MAX:
        return Verdict(BAD, f"ratio_above_thresh {ratio:.2f} too high", scores)
    if ratio < RATIO_ABOVE_THRESH_MIN:
        return Verdict(BAD, f"ratio_above_thresh {ratio:.2f} too low", scores)

    xavg = np.average(normalised, 0)
    yavg = np.average(normalised, 1)
    axis = list(range(HALF_DIMS))

    centre_x = np.average(axis, 0, xavg)
    centre_y = np.average(axis, 0, yavg)
    scores["centre_x"] = float(centre_x)
    scores["centre_y"] = float(centre_y)
    good_centre = not (
        (HALF_DIMS // 2 - centre_x > MAX_CENTER_SKEW)
        or (HALF_DIMS // 2 - centre_y > MAX_CENTER_SKEW)
    )

    high_x = np.argmax(np.cumsum(xavg) > (np.sum(xavg) * EDGE_THRESH))
    low_x = len(xavg) - np.argmax(np.cumsum(np.flip(xavg)) > (np.sum(xavg) * EDGE_THRESH))
    high_y = np.argmax(np.cumsum(yavg) > (np.sum(yavg) * EDGE_THRESH))
    low_y = len(yavg) - np.argmax(np.cumsum(np.flip(yavg)) > (np.sum(yavg) * EDGE_THRESH))
    # Not halved: the image is already downsampled by two in each dimension.
    radius = ((high_x - low_x) + (high_y - low_y)) / 2.0
    scores["radius"] = float(radius)
    good_radius = abs(radius - AVG_RADIUS) <= MAX_RADIUS_ERROR

    reference = ideal_disc_axis()
    gof_x = np.sum(np.abs(xavg - reference)) / HALF_DIMS
    gof_y = np.sum(np.abs(yavg - reference)) / HALF_DIMS
    scores["gof_x"] = float(gof_x)
    scores["gof_y"] = float(gof_y)
    good_fit = not (gof_x > MAX_GOF or gof_y > MAX_GOF)

    if good_centre and good_radius and good_fit:
        return Verdict(GOOD, None, scores)
    failed = [
        name
        for name, ok in (("centre", good_centre), ("radius", good_radius), ("fit", good_fit))
        if not ok
    ]
    return Verdict(BAD, "failed " + "+".join(failed), scores)


# -------------------------------------------------------------------- disc_v1

#: Angular samples in the polar resampling.  The measurements are all averages over
#: this axis, so more angles cost little and reduce noise.
DISC_ANGLES = 720
#: Working resolution.  Halving 1280 keeps the limb several pixels wide while making
#: the polar transform cheap.
DISC_SIZE = 640

#: Radial bands, in units of the expected solar radius, used to characterise the
#: profile.  The on-disc band avoids both the centre (where limb darkening and
#: filaments live) and the limb itself; the off-limb band sits outside the corona's
#: steepest falloff but inside the frame -- 1.45R of a ~193px radius is 280px, within
#: the 320px half-width.
DISC_ON_BAND = (0.30, 0.70)
DISC_OFF_BAND = (1.25, 1.45)
#: Where the limb is allowed to be found, again in units of expected radius.
DISC_SEARCH_BAND = (0.75, 1.25)


def disc_profile(image, expected_radius, size=DISC_SIZE, angles=DISC_ANGLES):
    """Measure the solar disc from its azimuthally averaged radial profile.

    Returns a dict of continuous measurements, or with None values where the frame
    is too degenerate to measure.  Separated from the verdict logic so that
    ``filter_FITS.py calibrate`` can gather the same numbers over a sample.

    Averaging over angle is the entire point.  An active region is localised in
    angle, so it barely shifts the averaged profile; the existing ``geometry_v1``
    averages along image rows and columns instead, which a bright region shifts
    directly -- the cause of its 36.9% rejection rate.

    `expected_radius` is in pixels *at the working resolution*, i.e. DIAM_SUN / 4.
    """
    empty = {"limb_contrast": None, "radius_ratio": None, "limb_width": None}
    if image is None or image.ndim != 2 or min(image.shape) < 16:
        return empty
    if not expected_radius or not np.isfinite(expected_radius) or expected_radius <= 0:
        return empty

    working = cv.resize(
        np.nan_to_num(np.asarray(image, dtype=np.float32), nan=0.0, posinf=0.0, neginf=0.0),
        (size, size),
        interpolation=cv.INTER_AREA,
    )
    # Radiance spans orders of magnitude; log compresses a flare so it cannot
    # dominate the profile the way it dominates a linear mean.
    working = np.log1p(np.clip(working, 0.0, None))

    max_radius = size / 2.0
    polar = cv.warpPolar(
        working,
        (int(max_radius), angles),
        (size / 2.0, size / 2.0),
        max_radius,
        cv.WARP_POLAR_LINEAR,
    )
    profile = polar.mean(axis=0)
    radii = np.arange(profile.size, dtype=np.float64)
    scale = expected_radius

    def band(low, high):
        mask = (radii > low * scale) & (radii < high * scale)
        return profile[mask] if mask.any() else np.array([])

    on_values, off_values = band(*DISC_ON_BAND), band(*DISC_OFF_BAND)
    if on_values.size < 3 or off_values.size < 3:
        return empty
    on_disc, off_limb = float(np.median(on_values)), float(np.median(off_values))
    if not np.isfinite(on_disc) or not np.isfinite(off_limb):
        return empty
    if on_disc <= 0:
        # No signal at all where the disc should be.  This is a conclusion, not a
        # failure to measure: report zero contrast so the caller rejects the frame
        # rather than abstaining on it.
        return {"limb_contrast": 0.0, "radius_ratio": None, "limb_width": None}

    contrast = (on_disc - off_limb) / on_disc
    measured = {"limb_contrast": float(contrast), "radius_ratio": None, "limb_width": None}
    if on_disc <= off_limb:
        # No radial falloff at all: there is no disc here to measure the size of.
        return measured

    # The limb is the steepest descent of the smoothed profile.  Taken on the
    # *averaged* profile rather than per angle: fitting each angle separately failed
    # outright on frames where too few angles yielded a usable edge.
    smoothed = cv.GaussianBlur(profile.astype(np.float32).reshape(1, -1), (9, 1), 0).ravel()
    derivative = np.gradient(smoothed.astype(np.float64))
    low, high = DISC_SEARCH_BAND
    window = (radii >= low * scale) & (radii <= high * scale)
    if window.sum() < 5:
        return measured
    indices = np.flatnonzero(window)
    local = derivative[indices]
    trough = int(indices[np.argmin(local)])
    depth = derivative[trough]
    if not np.isfinite(depth) or depth >= 0:
        return measured  # no descending edge: not a disc

    measured["radius_ratio"] = float(trough / scale)

    # Width of the descent at half its depth.  A disc displaced by d smears the
    # averaged limb across roughly 2d, so this is what makes the measurement
    # sensitive to decentring without ever having to locate a centre.
    half = depth / 2.0
    left = trough
    while left > 0 and derivative[left - 1] <= half:
        left -= 1
    right = trough
    last = derivative.size - 1
    while right < last and derivative[right + 1] <= half:
        right += 1
    measured["limb_width"] = float((right - left + 1) / scale)
    return measured


#: Per-band acceptance ranges for the disc_v1 measurements.
#:
#: Generated by ``filter_FITS.py calibrate``; regenerate and paste the output here
#: rather than hand-editing.  Derived from 400 frames per band sampled evenly across
#: the archive's whole time range, with 0/400 sample failures in every band.
#:
#: The bounds sit well outside the observed spread of good frames, because the
#: false-positive budget for this detector is <0.1%: the method it replaces was
#: discarding 36.9% of perfectly good frames, which is the entire reason it exists.
#: Only the diagnostic side of each measurement is bounded -- see CALIBRATION_SIDES
#: in filter_FITS.py.
#:
#: Note what the wide contrast floors mean in practice.  Good frames in most bands
#: reach as low as 0.11-0.22 contrast, so a floor tight enough to catch a partially
#: degraded frame would reject real data.  At this budget limb_contrast therefore
#: catches total signal loss (all-zero, NaN, uniform, which measure exactly 0.0) and
#: little else; size and shape faults are caught by the other two measurements.
#: 171A is the exception, where good frames never drop below 0.68.
#:
#: A measurement whose calibrated bounds cannot achieve the budget is set to None
#: here, which records it as a score without letting it fail a frame.
DISC_BOUNDS = {
    94: {
        "limb_contrast": (0.02, 1.0),
        "radius_ratio": (0.9631, 1.1176),
        "limb_width": (0.0, 0.1217),
    },
    131: {
        "limb_contrast": (0.02, 1.0),
        "radius_ratio": (0.9578, 1.0854),
        "limb_width": (0.0, 0.1011),
    },
    171: {
        "limb_contrast": (0.3183, 1.0),
        "radius_ratio": (0.9616, 1.1062),
        "limb_width": (0.0, 0.1375),
    },
    195: {
        "limb_contrast": (0.02, 1.0),
        "radius_ratio": (0.9676, 1.1167),
        "limb_width": (0.0, 0.1704),
    },
    284: {
        "limb_contrast": (0.02, 1.0),
        "radius_ratio": (0.9548, 1.1489),
        "limb_width": (0.0, 0.3177),
    },
    304: {
        "limb_contrast": (0.02, 1.0),
        "radius_ratio": (0.9548, 1.0813),
        "limb_width": (0.0, 0.1016),
    },
}


def disc_v1(image, wavelength, header=None, bounds=None):
    """Verify that the frame holds a correctly sized, centred, sharp solar disc.

    Replaces ``geometry_v1``'s intensity-weighted centroid, which a large active
    region drags off centre -- rejecting entire days of good data during exactly the
    solar activity most worth watching.

    Everything measured here is an average over angle, so where the bright regions
    sit on the disc does not move it.  What it deliberately cannot see is rotation:
    ``rotate`` and ``yaw_flip`` leave every azimuthally averaged quantity unchanged,
    and one frame carries no absolute rotation reference.  Those belong to
    ``temporal_v1`` and ``crosssat_v1``.

    The expected radius comes from the frame's own ``DIAM_SUN``, which varies by
    3.3% over the year with the Earth-Sun distance.
    """
    bounds = bounds or DISC_BOUNDS
    if image is None:
        return Verdict(BAD, "no image data", {})

    diameter = (header or {}).get("diam_sun")
    if not diameter:
        return Verdict(UNKNOWN, "no DIAM_SUN in header", {})

    # DIAM_SUN is a full diameter at full resolution; halve for radius, halve again
    # for the working downsample.
    expected_radius = float(diameter) / 4.0
    measured = disc_profile(image, expected_radius)
    scores = {k: v for k, v in measured.items() if v is not None}

    if measured["limb_contrast"] is None:
        return Verdict(UNKNOWN, "radial profile not measurable", scores)

    limits = bounds.get(wavelength)
    if limits is None:
        return Verdict(UNKNOWN, f"no disc bounds for band {wavelength}", scores)

    failed = []
    for name, limit in limits.items():
        if limit is None:
            continue  # recorded as a score, not used to judge
        value = measured.get(name)
        if value is None:
            failed.append(f"{name} unmeasurable")
        elif not (limit[0] <= value <= limit[1]):
            failed.append(f"{name} {value:.3f} outside [{limit[0]:.3f}, {limit[1]:.3f}]")

    if failed:
        return Verdict(BAD, "; ".join(failed), scores)
    return Verdict(GOOD, None, scores)


# ---------------------------------------------------------------- temporal_v1

#: Frames either side used to build the local reference.  At a 4-minute cadence,
#: +/-6 spans about 50 minutes -- long enough to be robust, short enough that real
#: solar evolution does not dominate.
TEMPORAL_WINDOW = 6
#: Robust z-score above which a frame's brightness is judged anomalous.
TEMPORAL_BRIGHTNESS_Z = 8.0
#: Structural difference from neighbours, relative to the neighbours' own churn.
TEMPORAL_STRUCTURE_RATIO = 6.0
#: Below this relative difference two consecutive frames are considered identical.
FROZEN_TOLERANCE = 1e-6


def _robust_z(values, index, window):
    """Median-absolute-deviation z-score of values[index] against its neighbours."""
    low = max(0, index - window)
    high = min(len(values), index + window + 1)
    neighbours = np.array(
        [values[i] for i in range(low, high) if i != index and np.isfinite(values[i])]
    )
    if neighbours.size < 3:
        return None
    median = np.median(neighbours)
    mad = np.median(np.abs(neighbours - median))
    if mad <= 0:
        # A perfectly steady neighbourhood: fall back to the spread, and treat an
        # exactly-constant one as uninformative rather than infinitely sensitive.
        spread = neighbours.std()
        if spread <= 0:
            return 0.0 if values[index] == median else np.inf
        return abs(values[index] - median) / spread
    return abs(values[index] - median) / (1.4826 * mad)


def temporal_v1(series, window=TEMPORAL_WINDOW, brightness_z=TEMPORAL_BRIGHTNESS_Z,
                structure_ratio=TEMPORAL_STRUCTURE_RATIO,
                frozen_tolerance=FROZEN_TOLERANCE):
    """Judge each frame against its neighbours in the same band and satellite.

    Catches what single-frame checks cannot: a frame that is individually plausible
    but inconsistent with the minutes either side of it -- a frozen duplicate, a
    brightness step, a substituted frame from another time.

    `series` must be ordered by time and come from one (satellite, wavelength)
    stream; gaps are fine, they simply widen the neighbourhood in wall-clock terms.
    Returns one Verdict per input frame.
    """
    count = len(series)
    means = [
        f.header.get("img_mean", np.nan) if not f.error else np.nan for f in series
    ]
    verdicts = []

    # Structural churn between consecutive thumbnails, when they are available.
    diffs = [np.nan] * count
    for i in range(1, count):
        a, b = series[i - 1].thumbnail, series[i].thumbnail
        if a is not None and b is not None and a.shape == b.shape:
            scale = float(np.abs(a).mean() + np.abs(b).mean()) / 2.0
            diffs[i] = float(np.abs(a - b).mean()) / scale if scale > 0 else np.nan

    for i, frame in enumerate(series):
        scores = {}
        if frame.error:
            verdicts.append(Verdict(BAD, f"unreadable: {frame.error}", {"read_ok": 0.0}))
            continue

        z = _robust_z(means, i, window)
        if z is not None:
            scores["brightness_z"] = float(z)

        # A frame identical to its predecessor means the feed stalled.
        if i > 0 and np.isfinite(diffs[i]) and diffs[i] < frozen_tolerance:
            scores["frame_diff"] = float(diffs[i])
            verdicts.append(Verdict(BAD, "identical to previous frame", scores))
            continue

        neighbour_diffs = np.array(
            [
                diffs[j]
                for j in range(max(1, i - window), min(count, i + window + 1))
                if j != i and np.isfinite(diffs[j])
            ]
        )
        if np.isfinite(diffs[i]) and neighbour_diffs.size >= 3:
            typical = float(np.median(neighbour_diffs))
            scores["frame_diff"] = float(diffs[i])
            scores["frame_diff_ratio"] = float(diffs[i] / typical) if typical > 0 else 0.0
            if typical > 0 and diffs[i] / typical > structure_ratio:
                verdicts.append(
                    Verdict(BAD, f"structural jump {diffs[i] / typical:.1f}x typical", scores)
                )
                continue

        if z is not None and z > brightness_z:
            verdicts.append(Verdict(BAD, f"brightness z={z:.1f}", scores))
            continue

        if z is None and not np.isfinite(diffs[i]):
            verdicts.append(Verdict(UNKNOWN, "insufficient context", scores))
            continue

        verdicts.append(Verdict(GOOD, None, scores))

    return verdicts


# ---------------------------------------------------------------- crosssat_v1

#: Relative disagreement between satellites above which one of them is wrong.
CROSSSAT_MAX_DIFF = 0.35
#: Alignment shift beyond which the two views cannot be meaningfully compared.
#: Geostationary parallax is under ~7px at 2.5 arcsec/px; more means a pointing fault.
CROSSSAT_MAX_SHIFT = 12.0
#: Plausible range for the calibration gain between the two flight models.  Wide
#: enough to absorb genuine instrument differences, narrow enough that a blackout or
#: a saturation blowout cannot be rescaled into agreement.
CROSSSAT_GAIN_RANGE = (0.2, 5.0)


def _gain_match(source, reference):
    """Least-squares gain and offset putting `source` on `reference`'s scale.

    The two instruments are different flight models with different responses, so a
    raw difference conflates calibration with genuine disagreement.

    Returns (matched, gain, offset).  The gain is returned rather than hidden
    because matching is otherwise *too* effective: a frame uniformly scaled by 1e-4
    fits perfectly after rescaling, so a detector that only looked at the residual
    would call a total blackout a match.  The size of the correction is itself
    evidence.
    """
    x = source.ravel().astype(np.float64)
    y = reference.ravel().astype(np.float64)
    finite = np.isfinite(x) & np.isfinite(y)
    if finite.sum() < 16:
        return source, 1.0, 0.0
    x, y = x[finite], y[finite]
    variance = float(((x - x.mean()) ** 2).sum())
    if variance <= 0:
        return source, 1.0, 0.0
    gain = float(((x - x.mean()) * (y - y.mean())).sum() / variance)
    offset = float(y.mean() - gain * x.mean())
    return source * gain + offset, gain, offset


def align_shift(a, b):
    """Sub-pixel translation between two thumbnails, via phase correlation."""
    if a is None or b is None or a.shape != b.shape:
        return None
    window = cv.createHanningWindow((a.shape[1], a.shape[0]), cv.CV_64F)
    (dx, dy), _ = cv.phaseCorrelate(
        np.nan_to_num(a.astype(np.float64)), np.nan_to_num(b.astype(np.float64)), window
    )
    return float(np.hypot(dx, dy))


def crosssat_v1(features_a, features_b, temporal_a=None, temporal_b=None,
                max_diff=CROSSSAT_MAX_DIFF, max_shift=CROSSSAT_MAX_SHIFT,
                gain_range=CROSSSAT_GAIN_RANGE):
    """Compare the two satellites' view of the same instant.

    GOES-16 and GOES-18 observe the same Sun from geostationary orbit, so at 1 AU
    their views differ only by a few pixels of parallax and by instrument
    calibration.  Sustained disagreement therefore means one of them is wrong.

    Deciding *which* one needs outside information: `temporal_a`/`temporal_b` are the
    corresponding :func:`temporal_v1` verdicts, and whichever frame also disagrees
    with its own history is blamed.  With no tiebreak available both are reported
    ``unknown`` rather than guessing -- a wrong attribution would discard a good
    frame and keep a bad one.

    Returns (verdict_a, verdict_b).
    """
    if features_a is None or features_b is None:
        return Verdict(UNKNOWN, "no counterpart", {}), Verdict(UNKNOWN, "no counterpart", {})
    if features_a.error or features_b.error:
        # Frame-level errors are not this detector's job; defer rather than double-count.
        return Verdict(UNKNOWN, "counterpart unreadable", {}), Verdict(
            UNKNOWN, "counterpart unreadable", {}
        )

    a, b = features_a.thumbnail, features_b.thumbnail
    if a is None or b is None or a.shape != b.shape:
        return Verdict(UNKNOWN, "no comparable thumbnail", {}), Verdict(
            UNKNOWN, "no comparable thumbnail", {}
        )

    scores = {}
    shift = align_shift(a, b)
    if shift is not None:
        scores["align_shift"] = shift

    matched, gain, offset = _gain_match(b, a)
    scale = float(np.abs(a).mean() + np.abs(matched).mean()) / 2.0
    difference = float(np.abs(a - matched).mean()) / scale if scale > 0 else np.inf
    scores["cross_diff"] = difference
    scores["cross_gain"] = gain
    scores["cross_offset"] = offset

    # A correction this large means the two frames are not on comparable scales at
    # all, whatever the residual says once it has been applied.
    gain_ok = gain_range[0] <= gain <= gain_range[1]
    agree = difference <= max_diff and gain_ok and (shift is None or shift <= max_shift)
    if agree:
        return Verdict(GOOD, None, scores), Verdict(GOOD, None, scores)

    reason = f"cross-satellite disagreement {difference:.2f}"
    if not gain_ok:
        reason = f"cross-satellite gain {gain:.3g} outside {gain_range}"
    elif shift is not None and shift > max_shift:
        reason = f"cross-satellite misalignment {shift:.1f}px"

    a_suspect = temporal_a is not None and temporal_a.verdict == BAD
    b_suspect = temporal_b is not None and temporal_b.verdict == BAD
    if a_suspect and not b_suspect:
        return Verdict(BAD, reason, scores), Verdict(GOOD, None, scores)
    if b_suspect and not a_suspect:
        return Verdict(GOOD, None, scores), Verdict(BAD, reason, scores)
    return Verdict(UNKNOWN, reason + " (source unclear)", scores), Verdict(
        UNKNOWN, reason + " (source unclear)", scores
    )


#: Detectors the drivers and bench expose by name.
FRAME_DETECTORS = {"header_v1": header_v1, "geometry_v1": geometry_v1, "disc_v1": disc_v1}
SERIES_DETECTORS = {"temporal_v1": temporal_v1, "crosssat_v1": crosssat_v1}
ALL_DETECTORS = tuple(FRAME_DETECTORS) + tuple(SERIES_DETECTORS)