remove-ai-watermarks/scripts/synthid_affine_lattice_probe.py

"""Probe a periodic carrier lattice and codeword with split confirmation.

The synchronization parameters are selected on one checkerboard of image
patches and scored on the disjoint checkerboard. The statistics measure complex
phase coherence and a content-whitened multichannel template match after
correcting every patch for its global origin. They do not use filenames,
metadata, or scene-class features.
"""

from __future__ import annotations

import json
import logging
import math
from dataclasses import asdict, dataclass, replace
from pathlib import Path
from typing import TYPE_CHECKING, Any

import click
import cv2
import numpy as np
from synthid_pixel_attack import jpeg_round_trip, load_rgb

from remove_ai_watermarks._synthid_confirmation import (
    registered_confirmation_passes as runtime_registered_confirmation_passes,
)
from remove_ai_watermarks._synthid_registered import RegisteredComponents, registered_components
from remove_ai_watermarks.synthid_detector import fold_residual_template, folded_template_score, unit_tile

if TYPE_CHECKING:
    from numpy.typing import NDArray

log = logging.getLogger(__name__)

FIXED_CANDIDATE_THRESHOLD = 0.28
PATCH_SHIFT_MIN_MARGIN = 0.45
PATCH_SHIFT_STRONG_MARGIN = 1.0
PATCH_SHIFT_MIN_Z = 2.5
OPPONENT_REGISTERED_FIXED_MIN = 0.16
OPPONENT_REGISTERED_RED_GREEN_MIN = 0.60
OPPONENT_REGISTERED_BLUE_YELLOW_MIN = 0.55
SAME_IMAGE_NULL_OFFSETS = (
    -2.0,
    -1.75,
    -1.5,
    -1.25,
    -1.0,
    -0.75,
    -0.5,
    -0.35,
    0.35,
    0.5,
    0.75,
    1.0,
    1.25,
    1.5,
    1.75,
    2.0,
)


def _webp_round_trip(pixels: NDArray[Any], quality: int) -> NDArray[Any]:
    """Apply one in-memory WebP encode/decode while returning RGB pixels."""
    if quality < 1 or quality > 101:
        raise ValueError("WebP quality must be between 1 and 101")
    success, encoded = cv2.imencode(
        ".webp",
        cv2.cvtColor(pixels, cv2.COLOR_RGB2BGR),
        [cv2.IMWRITE_WEBP_QUALITY, quality],
    )
    if not success:
        raise RuntimeError("WebP encoding failed")
    decoded = cv2.imdecode(encoded, cv2.IMREAD_COLOR)
    if decoded is None:
        raise RuntimeError("WebP decoding failed")
    return cv2.cvtColor(decoded, cv2.COLOR_BGR2RGB)


@dataclass(frozen=True)
class LatticeScore:
    """One select-confirm reciprocal-lattice observation."""

    selected_period: float
    selected_rotation_degrees: float
    selected_orientation_degrees: int
    selected_horizontal_reflection: bool
    selected_deskew_degrees: float
    deskew_candidate_count: int
    deskew_selection_coherence: float
    deskew_direct_selection_match: float
    deskew_direct_confirmation_match: float
    deskew_direct_joint_match: float
    peak_period: float
    peak_rotation_degrees: float
    peak_selection_coherence: float
    selection_coherence: float
    selection_candidate_p95: float
    selection_candidate_p99: float
    selection_excess_p95: float
    selection_excess_p99: float
    confirmation_coherence: float
    confirmation_candidate_p95: float
    confirmation_candidate_p99: float
    confirmation_excess_p95: float
    confirmation_excess_p99: float
    joint_coherence: float
    joint_excess_p99: float
    selected_shift_y: int
    selected_shift_x: int
    selection_codeword: float
    confirmation_codeword: float
    confirmation_codeword_shift_p95: float
    confirmation_codeword_shift_p99: float
    joint_codeword: float
    unknown_codeword_shift_y: int
    unknown_codeword_shift_x: int
    unknown_codeword_selection: float
    unknown_codeword_confirmation: float
    unknown_codeword_fixed_confirmation: float
    unknown_codeword_fixed_all: float
    unknown_codeword_confirmation_p95: float
    unknown_codeword_confirmation_p99: float
    unknown_codeword_excess_p99: float
    amplitude_shift_y: int
    amplitude_shift_x: int
    selection_amplitude: float
    confirmation_amplitude: float
    joint_amplitude: float
    selection_whitened_match: float
    confirmation_whitened_match: float
    joint_whitened_match: float
    amplitude_candidate_count: int
    amplitude_rerank_count: int
    canonical_template_score: float
    canonical_registered_template_score: float
    selection_patches: int
    confirmation_patches: int


@dataclass(frozen=True)
class _AmplitudeScore:
    """Amplitude-aware select-confirm score for one period candidate."""

    period: float
    selection: float
    confirmation: float
    shift_y: int
    shift_x: int
    unaligned_full: float
    aligned_full: float


@dataclass(frozen=True)
class SameImageNullScore:
    """Target-period coherence relative to neighboring periods in one image."""

    target_period: float
    target_selection_coherence: float
    target_confirmation_coherence: float
    joint_coherence: float
    selection_off_period_max: float
    confirmation_off_period_max: float
    joint_excess: float
    off_period_count: int


@dataclass(frozen=True)
class PatchShiftConsensusScore:
    """Robust per-patch agreement on one cyclic carrier phase."""

    period: float
    selected_shift_y: int
    selected_shift_x: int
    selection_trimmed_z: float
    confirmation_trimmed_z: float
    joint_trimmed_z: float
    selection_support_fraction: float
    confirmation_support_fraction: float
    joint_support_fraction: float
    selection_patches: int
    confirmation_patches: int


@dataclass(frozen=True)
class OpponentRegisteredScore:
    """Scale-registered opponent-color carrier observation."""

    selected_period: float
    spectral_period: float
    spectral_score: float
    fixed_score: float
    red_green_spatial: float
    blue_yellow_spatial: float
    candidate_count: int

    @property
    def decision_score(self) -> float:
        """Return the minimum normalized research gate margin."""
        return min(
            self.fixed_score / OPPONENT_REGISTERED_FIXED_MIN,
            self.red_green_spatial / OPPONENT_REGISTERED_RED_GREEN_MIN,
            self.blue_yellow_spatial / OPPONENT_REGISTERED_BLUE_YELLOW_MIN,
        )


def fixed_candidate_passes(score: float) -> bool:
    """Return the frozen precision-first fixed-period candidate verdict."""
    return math.isfinite(score) and score >= FIXED_CANDIDATE_THRESHOLD


def registered_confirmation_passes(score: LatticeScore) -> bool:
    """Return the frozen period-aware split-confirmation verdict."""
    return runtime_registered_confirmation_passes(
        score.selected_period,
        score.joint_coherence,
        score.joint_amplitude,
        score.unknown_codeword_fixed_confirmation,
    )


def patch_shift_recovery_passes(
    *,
    amplitude_margin: float,
    high_band_margin: float,
    periods_agree: bool,
    confirmation_passes: bool,
    joint_trimmed_z: float,
) -> bool:
    """Apply the frozen research-only content-adaptive recovery rule."""
    margins = (amplitude_margin, high_band_margin)
    return (
        all(math.isfinite(value) for value in (*margins, joint_trimmed_z))
        and periods_agree
        and confirmation_passes
        and min(margins) >= PATCH_SHIFT_MIN_MARGIN
        and max(margins) >= PATCH_SHIFT_STRONG_MARGIN
        and joint_trimmed_z >= PATCH_SHIFT_MIN_Z
    )


def _opponent_channels(values: NDArray[Any]) -> NDArray[Any]:
    """Return Green, Red-minus-Green, and Blue-minus-Yellow channels."""
    red = values[:, :, 0]
    green = values[:, :, 1]
    blue = values[:, :, 2]
    return np.stack((green, red - green, blue - 0.5 * (red + green)), axis=2)


def _opponent_color_pair(values: NDArray[Any]) -> NDArray[Any]:
    """Return the large-carrier Red-Green and Blue-Yellow planes."""
    red = values[:, :, 0]
    green = values[:, :, 1]
    blue = values[:, :, 2]
    return np.stack((red - green, blue - 0.5 * (red + green)), axis=2)


def template_harmonics(template: NDArray[Any], count: int = 16) -> tuple[NDArray[Any], NDArray[Any], NDArray[Any]]:
    """Return strong unique half-plane harmonics and their channel weights."""
    if template.ndim != 3 or template.shape[2] != 3:
        raise ValueError("template must have shape (height, width, 3)")
    if count < 1:
        raise ValueError("count must be positive")
    opponent = _opponent_channels(np.asarray(template, dtype=np.float64))
    spectrum = np.fft.fft2(opponent, axes=(0, 1))
    height, width = template.shape[:2]
    candidates: list[tuple[float, int, int, int, int]] = []
    for row in range(height):
        signed_row = row if row <= height // 2 else row - height
        for column in range(width):
            signed_column = column if column <= width // 2 else column - width
            if signed_row < 0 or (signed_row == 0 and signed_column <= 0):
                continue
            power = float(np.sum(np.abs(spectrum[row, column]) ** 2))
            candidates.append((power, signed_row, signed_column, row, column))
    candidates.sort(reverse=True)
    selected = candidates[:count]
    harmonics = np.asarray([(row, column) for _power, row, column, _y, _x in selected], dtype=np.float64)
    coefficients = np.asarray([spectrum[y, x] for _power, _row, _column, y, x in selected])
    weights = np.abs(coefficients)
    weight_sum = float(np.sum(weights))
    if weight_sum <= 0.0:
        raise ValueError("template has no nonzero periodic harmonics")
    coefficient_units = np.divide(
        coefficients,
        np.abs(coefficients),
        out=np.zeros_like(coefficients),
        where=np.abs(coefficients) > 1e-12,
    )
    return harmonics, weights / weight_sum, coefficient_units


def _patch_origins(height: int, width: int, patch_size: int, grid_size: int) -> list[tuple[int, int, int]]:
    if patch_size < 32 or height < patch_size or width < patch_size:
        raise ValueError("image must contain at least one patch of the requested size")
    if grid_size < 2:
        raise ValueError("grid size must be at least two")
    y_values = np.linspace(0, height - patch_size, min(grid_size, height // patch_size), dtype=np.int64)
    x_values = np.linspace(0, width - patch_size, min(grid_size, width // patch_size), dtype=np.int64)
    origins = []
    for y_index, y in enumerate(np.unique(y_values)):
        for x_index, x in enumerate(np.unique(x_values)):
            origins.append((int(y), int(x), (y_index + x_index) % 2))
    if {group for _y, _x, group in origins} != {0, 1}:
        raise ValueError("image geometry does not provide two patch groups")
    return origins


def _bilinear_sample(
    spectrum: NDArray[Any],
    y: NDArray[Any],
    x: NDArray[Any],
) -> NDArray[Any]:
    height, width = spectrum.shape
    y_floor = np.floor(y)
    x_floor = np.floor(x)
    y0 = y_floor.astype(np.int64) % height
    x0 = x_floor.astype(np.int64) % width
    y1 = (y0 + 1) % height
    x1 = (x0 + 1) % width
    dy = y - y_floor
    dx = x - x_floor
    return (
        spectrum[y0, x0] * (1.0 - dy) * (1.0 - dx)
        + spectrum[y1, x0] * dy * (1.0 - dx)
        + spectrum[y0, x1] * (1.0 - dy) * dx
        + spectrum[y1, x1] * dy * dx
    )


def _candidate_frequencies(
    periods: NDArray[Any],
    rotations_degrees: NDArray[Any],
    harmonics: NDArray[Any],
) -> tuple[NDArray[Any], NDArray[Any], NDArray[Any], NDArray[Any]]:
    period_grid, rotation_grid = np.meshgrid(periods, rotations_degrees, indexing="ij")
    flat_periods = period_grid.ravel()
    flat_rotations = rotation_grid.ravel()
    angles = np.deg2rad(flat_rotations)
    cosine = np.cos(angles)[:, None]
    sine = np.sin(angles)[:, None]
    base_y = harmonics[None, :, 0] / flat_periods[:, None]
    base_x = harmonics[None, :, 1] / flat_periods[:, None]
    frequencies_y = sine * base_x + cosine * base_y
    frequencies_x = cosine * base_x - sine * base_y
    return flat_periods, flat_rotations, frequencies_y, frequencies_x


def _patch_unit_values(
    pixels: NDArray[Any],
    origin_y: int,
    origin_x: int,
    patch_size: int,
    frequencies_y: NDArray[Any],
    frequencies_x: NDArray[Any],
) -> NDArray[Any]:
    patch = np.asarray(
        pixels[origin_y : origin_y + patch_size, origin_x : origin_x + patch_size],
        dtype=np.float32,
    )
    channels = _opponent_channels(patch)
    window_1d = np.hanning(patch_size).astype(np.float32)
    window = window_1d[:, None] * window_1d[None, :]
    sampled = np.empty((*frequencies_y.shape, 3), dtype=np.complex128)
    sample_y = frequencies_y * patch_size
    sample_x = frequencies_x * patch_size
    for channel in range(3):
        residual = channels[:, :, channel]
        residual = residual - cv2.GaussianBlur(
            residual,
            (0, 0),
            sigmaX=1.0,
            sigmaY=1.0,
            borderType=cv2.BORDER_REFLECT_101,
        )
        spectrum = np.fft.fft2(residual * window)
        sampled[:, :, channel] = _bilinear_sample(spectrum, sample_y, sample_x)
    phase_correction = np.exp(-2j * math.pi * (frequencies_y * origin_y + frequencies_x * origin_x))
    sampled *= phase_correction[:, :, None]
    magnitudes = np.abs(sampled)
    return np.divide(sampled, magnitudes, out=np.zeros_like(sampled), where=magnitudes > 1e-12)


def _coherence(unit_values: list[NDArray[Any]], weights: NDArray[Any]) -> NDArray[Any]:
    if not unit_values:
        raise ValueError("coherence needs at least one patch")
    values = np.stack(unit_values)
    coherence = np.abs(np.mean(values, axis=0))
    return np.sum(coherence * weights[None, :, :], axis=(1, 2))


def score_same_image_period_null(
    pixels: NDArray[Any],
    template: NDArray[Any],
    target_period: float,
    *,
    period_min: float = 7.5,
    period_max: float = 24.5,
    patch_size: int = 256,
    grid_size: int = 4,
    harmonic_count: int = 16,
) -> SameImageNullScore:
    """Compare a fixed target period with preregistered neighboring periods."""
    if not math.isfinite(target_period) or not period_min <= target_period <= period_max:
        raise ValueError("target period must be finite and inside the null range")
    off_periods = np.asarray(
        [
            target_period + offset
            for offset in SAME_IMAGE_NULL_OFFSETS
            if period_min <= target_period + offset <= period_max
        ],
        dtype=np.float64,
    )
    if not len(off_periods):
        raise ValueError("same-image period null needs at least one neighboring period")
    periods = np.concatenate((np.asarray([target_period]), off_periods))
    harmonics, weights, _coefficient_units = template_harmonics(template, harmonic_count)
    _periods, _rotations, frequencies_y, frequencies_x = _candidate_frequencies(
        periods,
        np.asarray([0.0]),
        harmonics,
    )
    grouped_values: dict[int, list[NDArray[Any]]] = {0: [], 1: []}
    for origin_y, origin_x, group in _patch_origins(*pixels.shape[:2], patch_size, grid_size):
        grouped_values[group].append(
            _patch_unit_values(
                pixels,
                origin_y,
                origin_x,
                patch_size,
                frequencies_y,
                frequencies_x,
            )
        )
    selection = _coherence(grouped_values[0], weights)
    confirmation = _coherence(grouped_values[1], weights)
    selection_off_max = float(np.max(selection[1:]))
    confirmation_off_max = float(np.max(confirmation[1:]))
    target_selection = float(selection[0])
    target_confirmation = float(confirmation[0])
    return SameImageNullScore(
        target_period=target_period,
        target_selection_coherence=target_selection,
        target_confirmation_coherence=target_confirmation,
        joint_coherence=min(target_selection, target_confirmation),
        selection_off_period_max=selection_off_max,
        confirmation_off_period_max=confirmation_off_max,
        joint_excess=min(
            target_selection - selection_off_max,
            target_confirmation - confirmation_off_max,
        ),
        off_period_count=len(off_periods),
    )


def _codeword_scores(
    selection_values: list[NDArray[Any]],
    confirmation_values: list[NDArray[Any]],
    selected_index: int,
    harmonics: NDArray[Any],
    coefficient_units: NDArray[Any],
    weights: NDArray[Any],
    tile_size: int,
) -> tuple[int, int, float, float, float, float]:
    selection_mean = np.mean(np.stack(selection_values), axis=0)[selected_index]
    confirmation_mean = np.mean(np.stack(confirmation_values), axis=0)[selected_index]
    shift_y, shift_x = np.meshgrid(np.arange(tile_size), np.arange(tile_size), indexing="ij")
    flat_y = shift_y.ravel()
    flat_x = shift_x.ravel()
    phase = np.exp(
        -2j * math.pi * (flat_y[:, None] * harmonics[None, :, 0] + flat_x[:, None] * harmonics[None, :, 1]) / tile_size
    )
    shifted_codewords = coefficient_units[None, :, :] * phase[:, :, None]

    def scores(values: NDArray[Any]) -> NDArray[Any]:
        agreement = np.real(values[None, :, :] * np.conj(shifted_codewords))
        return np.sum(agreement * weights[None, :, :], axis=(1, 2))

    selection_scores = scores(selection_mean)
    confirmation_scores = scores(confirmation_mean)
    selected_shift = int(np.argmax(selection_scores))
    return (
        int(flat_y[selected_shift]),
        int(flat_x[selected_shift]),
        float(selection_scores[selected_shift]),
        float(confirmation_scores[selected_shift]),
        float(np.quantile(confirmation_scores, 0.95)),
        float(np.quantile(confirmation_scores, 0.99)),
    )


def _unknown_codeword_scores(
    selection_values: list[NDArray[Any]],
    confirmation_values: list[NDArray[Any]],
    selected_index: int,
    harmonics: NDArray[Any],
    weights: NDArray[Any],
    tile_size: int,
) -> tuple[int, int, float, float, float, float, float, float]:
    """Fit an unknown patch codeword and confirm it on held-out harmonics."""
    if len(harmonics) < 4:
        raise ValueError("unknown-codeword confirmation needs at least four harmonics")
    selection_mean = np.mean(np.stack(selection_values), axis=0)[selected_index]
    confirmation_mean = np.mean(np.stack(confirmation_values), axis=0)[selected_index]
    cross_codeword = confirmation_mean * np.conj(selection_mean)

    shift_y, shift_x = np.meshgrid(np.arange(tile_size), np.arange(tile_size), indexing="ij")
    flat_y = shift_y.ravel()
    flat_x = shift_x.ravel()
    phase = np.exp(
        -2j * math.pi * (flat_y[:, None] * harmonics[None, :, 0] + flat_x[:, None] * harmonics[None, :, 1]) / tile_size
    )
    shifted_cross = cross_codeword[None, :, :] * phase[:, :, None]
    selection_mask = np.arange(len(harmonics)) % 2 == 0
    confirmation_mask = ~selection_mask

    def scores(mask: NDArray[Any]) -> NDArray[Any]:
        masked_weights = weights[mask]
        normalizer = float(np.sum(masked_weights))
        if normalizer <= 0.0:
            raise ValueError("unknown-codeword harmonic split has no template weight")
        return np.abs(np.sum(shifted_cross[:, mask, :] * masked_weights[None, :, :], axis=(1, 2))) / normalizer

    selection_scores = scores(selection_mask)
    confirmation_scores = scores(confirmation_mask)
    selected_shift = int(np.argmax(selection_scores))
    fixed_all = float(np.abs(np.sum(cross_codeword * weights)))
    return (
        int(flat_y[selected_shift]),
        int(flat_x[selected_shift]),
        float(selection_scores[selected_shift]),
        float(confirmation_scores[selected_shift]),
        float(confirmation_scores[0]),
        fixed_all,
        float(np.quantile(confirmation_scores, 0.95)),
        float(np.quantile(confirmation_scores, 0.99)),
    )


def _period_candidate_indices(
    periods: NDArray[Any],
    rotations_degrees: NDArray[Any],
    selection_scores: NDArray[Any],
    count: int,
) -> list[int]:
    """Return separated zero-rotation candidates ranked on selection patches."""
    candidates: list[int] = []
    eligible = np.flatnonzero(np.isclose(rotations_degrees, 0.0, atol=1e-12))
    for index in eligible[np.argsort(selection_scores[eligible])[::-1]]:
        if any(abs(float(periods[index] - periods[prior])) < 0.5 for prior in candidates):
            continue
        candidates.append(int(index))
        if len(candidates) == count:
            break
    if not candidates:
        raise ValueError("amplitude confirmation requires a zero-rotation candidate")
    return candidates


def _period_alias_candidate_indices(
    periods: NDArray[Any],
    rotations_degrees: NDArray[Any],
    base_indices: list[int],
) -> list[int]:
    """Expand coarse periods with octave aliases and adjacent grid bins."""
    eligible = np.flatnonzero(np.isclose(rotations_degrees, 0.0, atol=1e-12))
    eligible_periods = periods[eligible]
    candidates: list[int] = []
    seen: set[int] = set()
    for base_index in base_indices:
        for ratio in (1.0, 0.5, 2.0):
            target = float(periods[base_index] * ratio)
            if target < float(np.min(eligible_periods)) or target > float(np.max(eligible_periods)):
                continue
            center = int(np.argmin(np.abs(eligible_periods - target)))
            for position in range(max(0, center - 1), min(len(eligible), center + 2)):
                index = int(eligible[position])
                if index not in seen:
                    seen.add(index)
                    candidates.append(index)
    return candidates


def _cyclic_correlations(template: NDArray[Any], tile: NDArray[Any]) -> NDArray[Any]:
    template_spectrum = np.fft.fft2(template, axes=(0, 1))
    tile_spectrum = np.fft.fft2(tile, axes=(0, 1))
    return np.fft.ifft2(np.sum(template_spectrum * np.conj(tile_spectrum), axis=2)).real


def _canonical_pixels(pixels: NDArray[Any], template: NDArray[Any], period: float) -> NDArray[Any]:
    width = max(template.shape[1], round(pixels.shape[1] * template.shape[1] / period))
    height = max(template.shape[0], round(pixels.shape[0] * template.shape[0] / period))
    if (height, width) == pixels.shape[:2]:
        return pixels
    interpolation = cv2.INTER_AREA if width < pixels.shape[1] else cv2.INTER_CUBIC
    return np.asarray(cv2.resize(pixels, (width, height), interpolation=interpolation))


def _real_correlation(left: NDArray[Any], right: NDArray[Any]) -> float:
    denominator = float(np.linalg.norm(left) * np.linalg.norm(right))
    return float(np.real(np.vdot(right, left)) / denominator) if denominator > 0.0 else 0.0


def _opponent_period_curve(
    pixels: NDArray[Any],
    template: NDArray[Any],
    periods: NDArray[Any],
    *,
    harmonic_count: int,
) -> NDArray[Any]:
    """Return signed opponent-color template coherence over PERIODS."""
    template_opponent = _opponent_color_pair(np.asarray(template, dtype=np.float64))
    template_spectrum = np.fft.fft2(template_opponent, axes=(0, 1))
    power = np.sum(np.abs(template_spectrum) ** 2, axis=2)
    power[0, 0] = 0.0
    indices = np.argsort(power.ravel())[::-1][:harmonic_count]
    rows, columns = np.unravel_index(indices, power.shape)
    height, width = template.shape[:2]
    signed_rows = np.where(rows <= height // 2, rows, rows - height)
    signed_columns = np.where(columns <= width // 2, columns, columns - width)
    harmonics = np.column_stack((signed_rows, signed_columns)).astype(np.float64)
    coefficients = template_spectrum[rows, columns]

    image_height, image_width = pixels.shape[:2]
    sample_y = periods[:, None] ** -1 * harmonics[None, :, 0] * image_height
    sample_x = periods[:, None] ** -1 * harmonics[None, :, 1] * image_width
    sampled = np.empty((len(periods), len(harmonics), 2), dtype=np.complex128)
    image_opponent = _opponent_color_pair(np.asarray(pixels, dtype=np.float32))
    for channel in range(2):
        residual = image_opponent[:, :, channel]
        residual -= cv2.GaussianBlur(
            residual,
            (0, 0),
            sigmaX=1.0,
            sigmaY=1.0,
            borderType=cv2.BORDER_REFLECT_101,
        )
        sampled[:, :, channel] = _bilinear_sample(
            np.fft.fft2(residual),
            sample_y % image_height,
            sample_x % image_width,
        )
    numerator = np.real(np.sum(np.conj(coefficients)[None, :, :] * sampled, axis=(1, 2)))
    denominator = np.linalg.norm(coefficients) * np.linalg.norm(sampled, axis=(1, 2))
    return np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator > 0.0)


def score_opponent_registered(
    pixels: NDArray[Any],
    template: NDArray[Any],
    *,
    periods: NDArray[Any],
    harmonic_count: int = 30,
    candidate_count: int = 5,
    denoise_sigma: float = 1.0,
) -> OpponentRegisteredScore:
    """Search scale using the large-carrier opponent-color representation."""
    if pixels.ndim != 3 or pixels.shape[2] != 3:
        raise ValueError("pixels must have shape (height, width, 3)")
    if periods.ndim != 1 or not len(periods) or np.any(~np.isfinite(periods)) or np.any(periods <= 0.0):
        raise ValueError("periods must be a nonempty positive finite vector")
    if harmonic_count < 1 or candidate_count < 1:
        raise ValueError("harmonic and candidate counts must be positive")
    curve = _opponent_period_curve(
        pixels,
        template,
        periods,
        harmonic_count=harmonic_count,
    )
    rotations = np.zeros_like(periods)
    candidate_indices = _period_candidate_indices(periods, rotations, curve, candidate_count)
    observations: list[OpponentRegisteredScore] = []
    for index in candidate_indices:
        period = float(periods[index])
        canonical = _canonical_pixels(pixels, template, period)
        fixed_score, folded = folded_template_score(canonical, template, denoise_sigma)
        folded_opponent = _opponent_color_pair(folded)
        template_opponent = _opponent_color_pair(template)
        observations.append(
            OpponentRegisteredScore(
                selected_period=period,
                spectral_period=float(periods[int(np.argmax(curve))]),
                spectral_score=float(curve[index]),
                fixed_score=fixed_score,
                red_green_spatial=_real_correlation(folded_opponent[:, :, 0], template_opponent[:, :, 0]),
                blue_yellow_spatial=_real_correlation(folded_opponent[:, :, 1], template_opponent[:, :, 1]),
                candidate_count=len(candidate_indices),
            )
        )
    return max(observations, key=lambda score: score.decision_score)


def _content_whitened_patch_score(
    pixels: NDArray[Any],
    template: NDArray[Any],
    harmonics: NDArray[Any],
    *,
    denoise_sigma: float,
    noise_radius: int = 4,
    guard_radius: int = 1,
    ridge: float = 0.1,
) -> float:
    """Return a complex matched-filter cosine after local color whitening."""
    patch_height, patch_width = pixels.shape[:2]
    tile_height, tile_width = template.shape[:2]
    if patch_height % tile_height or patch_width % tile_width:
        raise ValueError("whitened patch dimensions must be multiples of the template dimensions")
    if noise_radius <= guard_radius or guard_radius < 0:
        raise ValueError("noise radius must exceed the nonnegative guard radius")
    if not math.isfinite(ridge) or ridge <= 0.0:
        raise ValueError("whitening ridge must be finite and positive")

    channels = _opponent_channels(np.asarray(pixels, dtype=np.float64))
    for channel in range(channels.shape[2]):
        residual = channels[:, :, channel]
        channels[:, :, channel] = residual - cv2.GaussianBlur(
            residual,
            (0, 0),
            sigmaX=denoise_sigma,
            sigmaY=denoise_sigma,
            borderType=cv2.BORDER_REFLECT_101,
        )
    spectrum = np.fft.fft2(channels, axes=(0, 1))
    template_spectrum = np.fft.fft2(_opponent_channels(np.asarray(template, dtype=np.float64)), axes=(0, 1))

    numerator = 0.0
    template_energy = 0.0
    observation_energy = 0.0
    for signed_row_value, signed_column_value in harmonics:
        signed_row = round(float(signed_row_value))
        signed_column = round(float(signed_column_value))
        bin_y = (signed_row * patch_height // tile_height) % patch_height
        bin_x = (signed_column * patch_width // tile_width) % patch_width
        noise = np.asarray(
            [
                spectrum[(bin_y + offset_y) % patch_height, (bin_x + offset_x) % patch_width]
                for offset_y in range(-noise_radius, noise_radius + 1)
                for offset_x in range(-noise_radius, noise_radius + 1)
                if max(abs(offset_y), abs(offset_x)) > guard_radius
            ]
        )
        covariance = noise.T @ np.conj(noise) / len(noise)
        local_power = float(np.trace(covariance).real / covariance.shape[0])
        covariance += np.eye(covariance.shape[0]) * (ridge * local_power + 1e-12)
        template_value = template_spectrum[signed_row % tile_height, signed_column % tile_width]
        observation = spectrum[bin_y, bin_x]
        whitened_template = np.linalg.solve(covariance, template_value)
        whitened_observation = np.linalg.solve(covariance, observation)
        numerator += float(np.vdot(template_value, whitened_observation).real)
        template_energy += float(np.vdot(template_value, whitened_template).real)
        observation_energy += float(np.vdot(observation, whitened_observation).real)
    denominator = math.sqrt(max(0.0, template_energy * observation_energy))
    return numerator / denominator if denominator > 1e-12 else 0.0


def _content_whitened_score(
    pixels: NDArray[Any],
    template: NDArray[Any],
    harmonics: NDArray[Any],
    *,
    period: float,
    patch_size: int,
    grid_size: int,
    denoise_sigma: float,
) -> tuple[float, float]:
    """Score a selected period on disjoint patch groups after local whitening."""
    canonical = _canonical_pixels(pixels, template, period)
    tile_height, tile_width = template.shape[:2]
    grouped_scores: dict[int, list[float]] = {0: [], 1: []}
    for origin_y, origin_x, group in _patch_origins(*canonical.shape[:2], patch_size, grid_size):
        aligned_y = (origin_y // tile_height) * tile_height
        aligned_x = (origin_x // tile_width) * tile_width
        patch = canonical[aligned_y : aligned_y + patch_size, aligned_x : aligned_x + patch_size]
        grouped_scores[group].append(
            _content_whitened_patch_score(
                patch,
                template,
                harmonics,
                denoise_sigma=denoise_sigma,
            )
        )
    return float(np.mean(grouped_scores[0])), float(np.mean(grouped_scores[1]))


def _amplitude_score(
    pixels: NDArray[Any],
    template: NDArray[Any],
    *,
    period: float,
    patch_size: int,
    grid_size: int,
    denoise_sigma: float,
) -> _AmplitudeScore:
    """Select cyclic phase on one patch group and confirm it on the other."""
    canonical = _canonical_pixels(pixels, template, period)
    grouped_units: dict[int, list[NDArray[Any]]] = {0: [], 1: []}
    tile_height, tile_width = template.shape[:2]
    for origin_y, origin_x, group in _patch_origins(*canonical.shape[:2], patch_size, grid_size):
        aligned_y = (origin_y // tile_height) * tile_height
        aligned_x = (origin_x // tile_width) * tile_width
        folded = fold_residual_template(
            canonical[aligned_y : aligned_y + patch_size, aligned_x : aligned_x + patch_size],
            tile_height=tile_height,
            tile_width=tile_width,
            denoise_sigma=denoise_sigma,
        )
        unit, _norm = unit_tile(folded)
        grouped_units[group].append(unit)
    selection_tile, _selection_norm = unit_tile(np.mean(grouped_units[0], axis=0))
    confirmation_tile, _confirmation_norm = unit_tile(np.mean(grouped_units[1], axis=0))
    selection_correlations = _cyclic_correlations(template, selection_tile)
    confirmation_correlations = _cyclic_correlations(template, confirmation_tile)
    selected_shift = int(np.argmax(selection_correlations))
    shift_y, shift_x = np.unravel_index(selected_shift, selection_correlations.shape)
    unaligned_full, full_folded = folded_template_score(canonical, template, denoise_sigma)
    aligned_full, _aligned_norm = unit_tile(np.roll(full_folded, shift=(shift_y, shift_x), axis=(0, 1)))
    return _AmplitudeScore(
        period=period,
        selection=float(selection_correlations[shift_y, shift_x]),
        confirmation=float(confirmation_correlations[shift_y, shift_x]),
        shift_y=int(shift_y),
        shift_x=int(shift_x),
        unaligned_full=unaligned_full,
        aligned_full=float(np.sum(template * aligned_full)),
    )


def _robust_shift_z(correlations: NDArray[Any]) -> NDArray[Any]:
    """Standardize each patch against its own cyclic-shift null."""
    flattened = correlations.reshape(correlations.shape[0], -1)
    center = np.median(flattened, axis=1, keepdims=True)
    mad = np.median(np.abs(flattened - center), axis=1, keepdims=True)
    scale = 1.4826 * mad
    fallback = np.std(flattened, axis=1, keepdims=True)
    scale = np.where(scale > 1e-12, scale, fallback)
    standardized = np.divide(
        flattened - center,
        scale,
        out=np.zeros_like(flattened),
        where=scale > 1e-12,
    )
    return standardized.reshape(correlations.shape)


def _top_half_mean(values: NDArray[Any], *, axis: int) -> NDArray[Any]:
    """Average the strongest half without letting one patch dominate."""
    count = values.shape[axis]
    retained = max(1, (count + 1) // 2)
    partitioned = np.partition(values, count - retained, axis=axis)
    indices = np.arange(count - retained, count)
    return np.mean(np.take(partitioned, indices, axis=axis), axis=axis)


def _shift_support_fraction(correlations: NDArray[Any], shift_y: int, shift_x: int) -> float:
    """Return the patch fraction placing one shift in its upper five percent."""
    selected = correlations[:, shift_y, shift_x]
    flattened = correlations.reshape(correlations.shape[0], -1)
    less = np.sum(flattened < selected[:, None], axis=1)
    equal = np.sum(flattened == selected[:, None], axis=1)
    percentile = (less + 0.5 * equal) / flattened.shape[1]
    return float(np.mean(percentile >= 0.95))


def score_patch_shift_consensus(
    pixels: NDArray[Any],
    template: NDArray[Any],
    period: float,
    *,
    patch_size: int = 256,
    grid_size: int = 4,
    denoise_sigma: float = 1.0,
) -> PatchShiftConsensusScore:
    """Select a robust cyclic phase and confirm it on disjoint patches.

    This probe targets a content-adaptive encoder that may place the shared
    detection carrier strongly in only part of an image. Every patch is
    normalized against its own 2-D cyclic-shift distribution before the
    strongest half are pooled. It is a research statistic, not a detector
    threshold.
    """
    if not math.isfinite(period) or period <= 0.0:
        raise ValueError("period must be finite and positive")
    canonical = _canonical_pixels(pixels, template, period)
    tile_height, tile_width = template.shape[:2]
    grouped_correlations: dict[int, list[NDArray[Any]]] = {0: [], 1: []}
    for origin_y, origin_x, group in _patch_origins(*canonical.shape[:2], patch_size, grid_size):
        aligned_y = (origin_y // tile_height) * tile_height
        aligned_x = (origin_x // tile_width) * tile_width
        folded = fold_residual_template(
            canonical[aligned_y : aligned_y + patch_size, aligned_x : aligned_x + patch_size],
            tile_height=tile_height,
            tile_width=tile_width,
            denoise_sigma=denoise_sigma,
        )
        unit, _norm = unit_tile(folded)
        grouped_correlations[group].append(_cyclic_correlations(template, unit))

    selection = np.stack(grouped_correlations[0])
    confirmation = np.stack(grouped_correlations[1])
    selection_z = _robust_shift_z(selection)
    confirmation_z = _robust_shift_z(confirmation)
    selection_trimmed = _top_half_mean(selection_z, axis=0)
    selected_shift = int(np.argmax(selection_trimmed))
    shift_y, shift_x = np.unravel_index(selected_shift, selection_trimmed.shape)
    selection_score = float(selection_trimmed[shift_y, shift_x])
    confirmation_score = float(_top_half_mean(confirmation_z[:, shift_y, shift_x], axis=0))
    selection_support = _shift_support_fraction(selection, int(shift_y), int(shift_x))
    confirmation_support = _shift_support_fraction(confirmation, int(shift_y), int(shift_x))
    return PatchShiftConsensusScore(
        period=period,
        selected_shift_y=int(shift_y),
        selected_shift_x=int(shift_x),
        selection_trimmed_z=selection_score,
        confirmation_trimmed_z=confirmation_score,
        joint_trimmed_z=min(selection_score, confirmation_score),
        selection_support_fraction=selection_support,
        confirmation_support_fraction=confirmation_support,
        joint_support_fraction=min(selection_support, confirmation_support),
        selection_patches=len(selection),
        confirmation_patches=len(confirmation),
    )


def score_lattice(
    pixels: NDArray[Any],
    template: NDArray[Any],
    *,
    periods: NDArray[Any],
    rotations_degrees: NDArray[Any],
    patch_size: int = 256,
    grid_size: int = 4,
    harmonic_count: int = 16,
    amplitude_candidate_count: int = 5,
    denoise_sigma: float = 1.0,
) -> LatticeScore:
    """Select a lattice on one patch group and confirm it on the other."""
    if pixels.ndim != 3 or pixels.shape[2] != 3:
        raise ValueError("pixels must have shape (height, width, 3)")
    if periods.ndim != 1 or not len(periods) or np.any(~np.isfinite(periods)) or np.any(periods <= 0.0):
        raise ValueError("periods must be a nonempty positive finite vector")
    if rotations_degrees.ndim != 1 or not len(rotations_degrees) or np.any(~np.isfinite(rotations_degrees)):
        raise ValueError("rotations must be a nonempty finite vector")
    if not math.isfinite(denoise_sigma) or denoise_sigma <= 0.0:
        raise ValueError("denoise sigma must be finite and positive")
    if amplitude_candidate_count < 1:
        raise ValueError("amplitude candidate count must be positive")
    harmonics, weights, coefficient_units = template_harmonics(template, harmonic_count)
    flat_periods, flat_rotations, frequencies_y, frequencies_x = _candidate_frequencies(
        periods,
        rotations_degrees,
        harmonics,
    )
    grouped_values: dict[int, list[NDArray[Any]]] = {0: [], 1: []}
    for origin_y, origin_x, group in _patch_origins(*pixels.shape[:2], patch_size, grid_size):
        grouped_values[group].append(
            _patch_unit_values(
                pixels,
                origin_y,
                origin_x,
                patch_size,
                frequencies_y,
                frequencies_x,
            )
        )
    selection_scores = _coherence(grouped_values[0], weights)
    confirmation_scores = _coherence(grouped_values[1], weights)
    peak_index = int(np.argmax(selection_scores))
    coarse_candidate_indices = _period_candidate_indices(
        flat_periods,
        flat_rotations,
        selection_scores,
        amplitude_candidate_count,
    )
    candidate_indices = _period_alias_candidate_indices(
        flat_periods,
        flat_rotations,
        coarse_candidate_indices,
    )
    whitened_scores = [
        _content_whitened_score(
            pixels,
            template,
            harmonics,
            period=float(flat_periods[index]),
            patch_size=patch_size,
            grid_size=grid_size,
            denoise_sigma=denoise_sigma,
        )
        for index in candidate_indices
    ]
    whitened_order = np.argsort([selection for selection, _confirmation in whitened_scores])[::-1]
    rerank_candidates = [int(index) for index in whitened_order[:3]]
    rerank_amplitudes = [
        _amplitude_score(
            pixels,
            template,
            period=float(flat_periods[candidate_indices[index]]),
            patch_size=patch_size,
            grid_size=grid_size,
            denoise_sigma=denoise_sigma,
        )
        for index in rerank_candidates
    ]
    rerank_scores = [
        min(*whitened_scores[index], amplitude.selection, amplitude.confirmation)
        for index, amplitude in zip(rerank_candidates, rerank_amplitudes, strict=True)
    ]
    rerank_winner = int(np.argmax(rerank_scores))
    selected_candidate = rerank_candidates[rerank_winner]
    selected_index = candidate_indices[selected_candidate]
    amplitude = rerank_amplitudes[rerank_winner]
    selection_whitened, confirmation_whitened = whitened_scores[selected_candidate]
    selection_p95 = float(np.quantile(selection_scores, 0.95))
    selection_p99 = float(np.quantile(selection_scores, 0.99))
    confirmation_p95 = float(np.quantile(confirmation_scores, 0.95))
    confirmation_p99 = float(np.quantile(confirmation_scores, 0.99))
    confirmation = float(confirmation_scores[selected_index])
    selection = float(selection_scores[selected_index])
    shift_y, shift_x, selection_codeword, confirmation_codeword, codeword_p95, codeword_p99 = _codeword_scores(
        grouped_values[0],
        grouped_values[1],
        selected_index,
        harmonics,
        coefficient_units,
        weights,
        template.shape[0],
    )
    (
        unknown_shift_y,
        unknown_shift_x,
        unknown_selection,
        unknown_confirmation,
        unknown_fixed_confirmation,
        unknown_fixed_all,
        unknown_p95,
        unknown_p99,
    ) = _unknown_codeword_scores(
        grouped_values[0],
        grouped_values[1],
        selected_index,
        harmonics,
        weights,
        template.shape[0],
    )
    return LatticeScore(
        selected_period=float(flat_periods[selected_index]),
        selected_rotation_degrees=float(flat_rotations[selected_index]),
        selected_orientation_degrees=0,
        selected_horizontal_reflection=False,
        selected_deskew_degrees=0.0,
        deskew_candidate_count=0,
        deskew_selection_coherence=0.0,
        deskew_direct_selection_match=0.0,
        deskew_direct_confirmation_match=0.0,
        deskew_direct_joint_match=0.0,
        peak_period=float(flat_periods[peak_index]),
        peak_rotation_degrees=float(flat_rotations[peak_index]),
        peak_selection_coherence=float(selection_scores[peak_index]),
        selection_coherence=selection,
        selection_candidate_p95=selection_p95,
        selection_candidate_p99=selection_p99,
        selection_excess_p95=selection - selection_p95,
        selection_excess_p99=selection - selection_p99,
        confirmation_coherence=confirmation,
        confirmation_candidate_p95=confirmation_p95,
        confirmation_candidate_p99=confirmation_p99,
        confirmation_excess_p95=confirmation - confirmation_p95,
        confirmation_excess_p99=confirmation - confirmation_p99,
        joint_coherence=min(selection, confirmation),
        joint_excess_p99=min(selection - selection_p99, confirmation - confirmation_p99),
        selected_shift_y=shift_y,
        selected_shift_x=shift_x,
        selection_codeword=selection_codeword,
        confirmation_codeword=confirmation_codeword,
        confirmation_codeword_shift_p95=codeword_p95,
        confirmation_codeword_shift_p99=codeword_p99,
        joint_codeword=min(selection_codeword, confirmation_codeword),
        unknown_codeword_shift_y=unknown_shift_y,
        unknown_codeword_shift_x=unknown_shift_x,
        unknown_codeword_selection=unknown_selection,
        unknown_codeword_confirmation=unknown_confirmation,
        unknown_codeword_fixed_confirmation=unknown_fixed_confirmation,
        unknown_codeword_fixed_all=unknown_fixed_all,
        unknown_codeword_confirmation_p95=unknown_p95,
        unknown_codeword_confirmation_p99=unknown_p99,
        unknown_codeword_excess_p99=unknown_confirmation - unknown_p99,
        amplitude_shift_y=amplitude.shift_y,
        amplitude_shift_x=amplitude.shift_x,
        selection_amplitude=amplitude.selection,
        confirmation_amplitude=amplitude.confirmation,
        joint_amplitude=min(amplitude.selection, amplitude.confirmation),
        selection_whitened_match=selection_whitened,
        confirmation_whitened_match=confirmation_whitened,
        joint_whitened_match=min(selection_whitened, confirmation_whitened),
        amplitude_candidate_count=len(candidate_indices),
        amplitude_rerank_count=len(rerank_candidates),
        canonical_template_score=amplitude.unaligned_full,
        canonical_registered_template_score=amplitude.aligned_full,
        selection_patches=len(grouped_values[0]),
        confirmation_patches=len(grouped_values[1]),
    )


def score_orientation_bank(
    pixels: NDArray[Any],
    template: NDArray[Any],
    *,
    periods: NDArray[Any],
    rotations_degrees: NDArray[Any],
    patch_size: int = 256,
    grid_size: int = 4,
    harmonic_count: int = 16,
    amplitude_candidate_count: int = 5,
    denoise_sigma: float = 1.0,
) -> LatticeScore:
    """Select a right-angle orientation on selection-patch whitened match."""
    candidates = [
        score_lattice(
            np.rot90(pixels, k=quarter_turns),
            template,
            periods=periods,
            rotations_degrees=rotations_degrees,
            patch_size=patch_size,
            grid_size=grid_size,
            harmonic_count=harmonic_count,
            amplitude_candidate_count=amplitude_candidate_count,
            denoise_sigma=denoise_sigma,
        )
        for quarter_turns in range(4)
    ]
    selected_index = int(np.argmax([candidate.selection_whitened_match for candidate in candidates]))
    return replace(candidates[selected_index], selected_orientation_degrees=selected_index * 90)


def score_dihedral_bank(
    pixels: NDArray[Any],
    template: NDArray[Any],
    *,
    periods: NDArray[Any],
    rotations_degrees: NDArray[Any],
    patch_size: int = 256,
    grid_size: int = 4,
    harmonic_count: int = 16,
    amplitude_candidate_count: int = 5,
    denoise_sigma: float = 1.0,
) -> LatticeScore:
    """Select a rotation and optional reflection on selection patches."""
    candidates: list[LatticeScore] = []
    transforms: list[tuple[int, bool]] = []
    for reflected in (False, True):
        reflected_pixels = np.fliplr(pixels) if reflected else pixels
        for quarter_turns in range(4):
            candidates.append(
                score_lattice(
                    np.rot90(reflected_pixels, k=quarter_turns),
                    template,
                    periods=periods,
                    rotations_degrees=rotations_degrees,
                    patch_size=patch_size,
                    grid_size=grid_size,
                    harmonic_count=harmonic_count,
                    amplitude_candidate_count=amplitude_candidate_count,
                    denoise_sigma=denoise_sigma,
                )
            )
            transforms.append((quarter_turns * 90, reflected))
    selected_index = int(np.argmax([candidate.selection_whitened_match for candidate in candidates]))
    selected_orientation, selected_reflection = transforms[selected_index]
    return replace(
        candidates[selected_index],
        selected_orientation_degrees=selected_orientation,
        selected_horizontal_reflection=selected_reflection,
    )


def _rotate_fixed_canvas(pixels: NDArray[Any], angle_degrees: float) -> NDArray[Any]:
    """Rotate RGB pixels around the image center without changing the canvas."""
    height, width = pixels.shape[:2]
    transform = cv2.getRotationMatrix2D(((width - 1) / 2.0, (height - 1) / 2.0), angle_degrees, 1.0)
    return np.asarray(
        cv2.warpAffine(
            pixels,
            transform,
            (width, height),
            flags=cv2.INTER_CUBIC,
            borderMode=cv2.BORDER_REFLECT_101,
        )
    )


def _select_deskew_angle(
    pixels: NDArray[Any],
    template: NDArray[Any],
    periods: NDArray[Any],
    deskew_degrees: NDArray[Any],
    *,
    patch_size: int,
    grid_size: int,
    harmonic_count: int,
) -> tuple[float, float, float]:
    """Select one deskew angle by carrier coherence on selection patches."""
    harmonics, weights, _coefficient_units = template_harmonics(template, harmonic_count)
    flat_periods, flat_rotations, frequencies_y, frequencies_x = _candidate_frequencies(
        periods,
        deskew_degrees,
        harmonics,
    )
    selection_values = [
        _patch_unit_values(
            pixels,
            origin_y,
            origin_x,
            patch_size,
            frequencies_y,
            frequencies_x,
        )
        for origin_y, origin_x, group in _patch_origins(*pixels.shape[:2], patch_size, grid_size)
        if group == 0
    ]
    selection_scores = _coherence(selection_values, weights)
    selected_index = int(np.argmax(selection_scores))
    return (
        float(flat_periods[selected_index]),
        float(flat_rotations[selected_index]),
        float(selection_scores[selected_index]),
    )


def _affine_whitened_patch_score(
    pixels: NDArray[Any],
    expected_template: NDArray[Any],
    frequencies_y: NDArray[Any],
    frequencies_x: NDArray[Any],
    *,
    origin_y: int,
    origin_x: int,
    denoise_sigma: float,
    noise_radius: int = 4,
    guard_radius: int = 1,
    ridge: float = 0.1,
) -> float:
    """Match an affine carrier directly, without resampling the image."""
    patch_size = pixels.shape[0]
    channels = _opponent_channels(np.asarray(pixels, dtype=np.float64))
    for channel in range(channels.shape[2]):
        residual = channels[:, :, channel]
        channels[:, :, channel] = residual - cv2.GaussianBlur(
            residual,
            (0, 0),
            sigmaX=denoise_sigma,
            sigmaY=denoise_sigma,
            borderType=cv2.BORDER_REFLECT_101,
        )
    window_1d = np.hanning(patch_size)
    spectrum = np.fft.fft2(channels * (window_1d[:, None] * window_1d[None, :])[:, :, None], axes=(0, 1))
    offset_y, offset_x = zip(
        *(
            (y, x)
            for y in range(-noise_radius, noise_radius + 1)
            for x in range(-noise_radius, noise_radius + 1)
            if max(abs(y), abs(x)) > guard_radius
        ),
        strict=True,
    )
    offset_y_values = np.asarray(offset_y, dtype=np.float64)
    offset_x_values = np.asarray(offset_x, dtype=np.float64)
    sample_y = frequencies_y * patch_size
    sample_x = frequencies_x * patch_size
    phase_correction = np.exp(-2j * math.pi * (frequencies_y * origin_y + frequencies_x * origin_x))
    numerator = 0.0
    template_energy = 0.0
    observation_energy = 0.0
    for harmonic_index in range(len(frequencies_y)):
        y = np.asarray([sample_y[harmonic_index]])
        x = np.asarray([sample_x[harmonic_index]])
        observation = np.asarray(
            [_bilinear_sample(spectrum[:, :, channel], y, x)[0] for channel in range(channels.shape[2])]
        )
        observation *= phase_correction[harmonic_index]
        noise = np.stack(
            [
                _bilinear_sample(
                    spectrum[:, :, channel],
                    y + offset_y_values,
                    x + offset_x_values,
                )
                for channel in range(channels.shape[2])
            ],
            axis=1,
        )
        covariance = noise.T @ np.conj(noise) / len(noise)
        local_power = float(np.trace(covariance).real / covariance.shape[0])
        covariance += np.eye(covariance.shape[0]) * (ridge * local_power + 1e-12)
        template_value = expected_template[harmonic_index]
        whitened_template = np.linalg.solve(covariance, template_value)
        whitened_observation = np.linalg.solve(covariance, observation)
        numerator += float(np.vdot(template_value, whitened_observation).real)
        template_energy += float(np.vdot(template_value, whitened_template).real)
        observation_energy += float(np.vdot(observation, whitened_observation).real)
    denominator = math.sqrt(max(0.0, template_energy * observation_energy))
    return numerator / denominator if denominator > 1e-12 else 0.0


def _affine_whitened_score(
    pixels: NDArray[Any],
    template: NDArray[Any],
    harmonics: NDArray[Any],
    *,
    period: float,
    deskew_degrees: float,
    patch_size: int,
    grid_size: int,
    denoise_sigma: float,
) -> tuple[float, float]:
    """Score a rotated carrier in the original pixels on disjoint patches."""
    _periods, _rotations, frequencies_y, frequencies_x = _candidate_frequencies(
        np.asarray([period]),
        np.asarray([deskew_degrees]),
        harmonics,
    )
    height, width = pixels.shape[:2]
    input_rotation = -deskew_degrees
    forward = cv2.getRotationMatrix2D(
        ((width - 1) / 2.0, (height - 1) / 2.0),
        input_rotation,
        1.0,
    )
    inverse = cv2.invertAffineTransform(forward)
    origin_x, origin_y = inverse[:, 2]
    phase = np.exp(2j * math.pi * (harmonics[:, 0] * origin_y / period + harmonics[:, 1] * origin_x / period))
    template_spectrum = np.fft.fft2(_opponent_channels(np.asarray(template, dtype=np.float64)), axes=(0, 1))
    expected_template = (
        np.asarray(
            [
                template_spectrum[round(float(row)) % template.shape[0], round(float(column)) % template.shape[1]]
                for row, column in harmonics
            ]
        )
        * phase[:, None]
    )
    grouped_scores: dict[int, list[float]] = {0: [], 1: []}
    for patch_y, patch_x, group in _patch_origins(height, width, patch_size, grid_size):
        patch = pixels[patch_y : patch_y + patch_size, patch_x : patch_x + patch_size]
        grouped_scores[group].append(
            _affine_whitened_patch_score(
                patch,
                expected_template,
                frequencies_y[0],
                frequencies_x[0],
                origin_y=patch_y,
                origin_x=patch_x,
                denoise_sigma=denoise_sigma,
            )
        )
    return float(np.mean(grouped_scores[0])), float(np.mean(grouped_scores[1]))


def score_deskew_bank(
    pixels: NDArray[Any],
    template: NDArray[Any],
    *,
    periods: NDArray[Any],
    deskew_degrees: NDArray[Any],
    patch_size: int = 256,
    grid_size: int = 4,
    harmonic_count: int = 16,
    amplitude_candidate_count: int = 5,
    denoise_sigma: float = 1.0,
) -> LatticeScore:
    """Select a small-angle deskew operation on selection patches."""
    if deskew_degrees.ndim != 1 or not len(deskew_degrees) or np.any(~np.isfinite(deskew_degrees)):
        raise ValueError("deskew angles must be a nonempty finite vector")
    harmonics, _weights, _coefficient_units = template_harmonics(template, harmonic_count)
    selected_period, selected_angle, selection_coherence = _select_deskew_angle(
        pixels,
        template,
        periods,
        deskew_degrees,
        patch_size=patch_size,
        grid_size=grid_size,
        harmonic_count=harmonic_count,
    )
    direct_selection, direct_confirmation = _affine_whitened_score(
        pixels,
        template,
        harmonics,
        period=selected_period,
        deskew_degrees=selected_angle,
        patch_size=patch_size,
        grid_size=grid_size,
        denoise_sigma=denoise_sigma,
    )
    score = score_lattice(
        _rotate_fixed_canvas(pixels, selected_angle),
        template,
        periods=periods,
        rotations_degrees=np.asarray([0.0]),
        patch_size=patch_size,
        grid_size=grid_size,
        harmonic_count=harmonic_count,
        amplitude_candidate_count=amplitude_candidate_count,
        denoise_sigma=denoise_sigma,
    )
    return replace(
        score,
        selected_deskew_degrees=selected_angle,
        deskew_candidate_count=len(periods) * len(deskew_degrees),
        deskew_selection_coherence=selection_coherence,
        deskew_direct_selection_match=direct_selection,
        deskew_direct_confirmation_match=direct_confirmation,
        deskew_direct_joint_match=min(direct_selection, direct_confirmation),
    )


def _load_template(path: Path) -> NDArray[Any]:
    with np.load(path, allow_pickle=False) as artifact:
        template = np.asarray(artifact["template"], dtype=np.float64)
    if template.shape != (16, 16, 3) or not np.all(np.isfinite(template)):
        raise ValueError("template artifact does not contain a finite 16x16 RGB template")
    return template


@click.command()
@click.argument("template_path", type=click.Path(exists=True, dir_okay=False, path_type=Path))
@click.argument("images", nargs=-1, required=True, type=click.Path(exists=True, dir_okay=False, path_type=Path))
@click.option("--period-min", type=click.FloatRange(min=1.0), default=7.5, show_default=True)
@click.option("--period-max", type=click.FloatRange(min=1.0), default=24.5, show_default=True)
@click.option("--period-step", type=click.FloatRange(min=0.01), default=0.1, show_default=True)
@click.option("--rotation-min", type=float, default=-3.0, show_default=True)
@click.option("--rotation-max", type=float, default=3.0, show_default=True)
@click.option("--rotation-step", type=click.FloatRange(min=0.01), default=0.25, show_default=True)
@click.option("--input-angle", type=float, default=0.0, show_default=True)
@click.option("--deskew-search/--no-deskew-search", default=False, show_default=True)
@click.option("--input-scale", type=click.FloatRange(min=0.01), default=1.0, show_default=True)
@click.option("--crop-fraction", type=click.FloatRange(min=0.0, max=0.49), default=0.0, show_default=True)
@click.option("--jpeg-quality", type=click.IntRange(min=1, max=100), default=None)
@click.option("--webp-quality", type=click.IntRange(min=1, max=101), default=None)
@click.option("--input-rotation", type=click.Choice(("0", "90", "180", "270")), default="0", show_default=True)
@click.option(
    "--input-flip",
    type=click.Choice(("none", "horizontal", "vertical", "both")),
    default="none",
    show_default=True,
)
@click.option("--orientation-search/--no-orientation-search", default=False, show_default=True)
@click.option("--dihedral-search/--no-dihedral-search", default=False, show_default=True)
@click.option("--registered-period/--no-registered-period", default=False, show_default=True)
@click.option("--same-image-null/--no-same-image-null", default=False, show_default=True)
@click.option("--patch-shift-consensus/--no-patch-shift-consensus", default=False, show_default=True)
@click.option("--opponent-registered/--no-opponent-registered", default=False, show_default=True)
@click.option("--amplitude-candidate-count", type=click.IntRange(min=1), default=5, show_default=True)
@click.option("--report-out", type=click.Path(dir_okay=False, path_type=Path), required=True)
def main(
    template_path: Path,
    images: tuple[Path, ...],
    period_min: float,
    period_max: float,
    period_step: float,
    rotation_min: float,
    rotation_max: float,
    rotation_step: float,
    input_angle: float,
    deskew_search: bool,
    input_scale: float,
    crop_fraction: float,
    jpeg_quality: int | None,
    webp_quality: int | None,
    input_rotation: str,
    input_flip: str,
    orientation_search: bool,
    dihedral_search: bool,
    registered_period: bool,
    same_image_null: bool,
    patch_shift_consensus: bool,
    opponent_registered: bool,
    amplitude_candidate_count: int,
    report_out: Path,
) -> None:
    """Score IMAGES with split-confirmed reciprocal-lattice coherence."""
    logging.basicConfig(level=logging.INFO, format="%(message)s")
    if period_max < period_min or rotation_max < rotation_min:
        raise click.BadParameter("search maximum must be at least its minimum")
    if jpeg_quality is not None and webp_quality is not None:
        raise click.UsageError("JPEG and WebP transforms are mutually exclusive")
    enabled_geometric_searches = sum((orientation_search, dihedral_search, deskew_search))
    if enabled_geometric_searches > 1:
        raise click.UsageError("orientation, dihedral, and deskew search are mutually exclusive")
    if registered_period and enabled_geometric_searches:
        raise click.UsageError("registered-period scoring cannot be combined with a geometric search")
    if same_image_null and not registered_period:
        raise click.UsageError("same-image-null scoring requires registered-period scoring")
    if patch_shift_consensus and not registered_period:
        raise click.UsageError("patch-shift consensus requires registered-period scoring")
    periods = np.arange(period_min, period_max + period_step / 2.0, period_step)
    rotations = np.arange(rotation_min, rotation_max + rotation_step / 2.0, rotation_step)
    template = _load_template(template_path)
    rows = []
    for path in images:
        registered: RegisteredComponents | None = None
        period_null: SameImageNullScore | None = None
        patch_consensus: PatchShiftConsensusScore | None = None
        opponent_score: OpponentRegisteredScore | None = None
        try:
            pixels = load_rgb(path)
            if input_scale != 1.0:
                scaled_width = max(1, round(pixels.shape[1] * input_scale))
                scaled_height = max(1, round(pixels.shape[0] * input_scale))
                interpolation = cv2.INTER_AREA if input_scale < 1.0 else cv2.INTER_CUBIC
                pixels = cv2.resize(pixels, (scaled_width, scaled_height), interpolation=interpolation)
            if crop_fraction:
                crop_y = round(pixels.shape[0] * crop_fraction)
                crop_x = round(pixels.shape[1] * crop_fraction)
                pixels = pixels[crop_y:, crop_x:]
            if jpeg_quality is not None:
                pixels = jpeg_round_trip(pixels, jpeg_quality)
            if webp_quality is not None:
                pixels = _webp_round_trip(pixels, webp_quality)
            if input_angle:
                pixels = _rotate_fixed_canvas(pixels, input_angle)
            quarter_turns = int(input_rotation) // 90
            if quarter_turns:
                pixels = np.rot90(pixels, k=-quarter_turns)
            if input_flip in {"horizontal", "both"}:
                pixels = np.fliplr(pixels)
            if input_flip in {"vertical", "both"}:
                pixels = np.flipud(pixels)
            image_periods = periods
            image_rotations = rotations
            if registered_period:
                registered = registered_components(pixels, template, 1.0)
                image_periods = np.asarray([registered.selected_period])
                image_rotations = np.asarray([0.0])
                if same_image_null:
                    period_null = score_same_image_period_null(
                        pixels,
                        template,
                        registered.selected_period,
                    )
                if patch_shift_consensus:
                    patch_consensus = score_patch_shift_consensus(
                        pixels,
                        template,
                        registered.selected_period,
                    )
            if opponent_registered:
                opponent_score = score_opponent_registered(
                    pixels,
                    template,
                    periods=periods,
                )
            if dihedral_search:
                scorer = score_dihedral_bank
            elif orientation_search:
                scorer = score_orientation_bank
            elif deskew_search:
                score = score_deskew_bank(
                    pixels,
                    template,
                    periods=periods,
                    deskew_degrees=rotations,
                    amplitude_candidate_count=amplitude_candidate_count,
                )
                scorer = None
            else:
                scorer = score_lattice
            if scorer is not None:
                score = scorer(
                    pixels,
                    template,
                    periods=image_periods,
                    rotations_degrees=image_rotations,
                    amplitude_candidate_count=amplitude_candidate_count,
                )
        except (OSError, ValueError) as error:
            rows.append({"path": str(path), "status": "unsupported", "reason": str(error)})
            log.warning("%s: unsupported: %s", path, error)
            continue
        row: dict[str, Any] = {"path": str(path), "status": "scored", "score": asdict(score)}
        if registered is not None:
            row["registered"] = {**asdict(registered), "decision_score": registered.decision_score}
            row["registered_confirmation_passes"] = registered_confirmation_passes(score)
        if period_null is not None:
            row["same_image_null"] = asdict(period_null)
        if patch_consensus is not None:
            row["patch_shift_consensus"] = asdict(patch_consensus)
        if opponent_score is not None:
            row["opponent_registered"] = {
                **asdict(opponent_score),
                "decision_score": opponent_score.decision_score,
            }
        rows.append(row)
        log.info(
            "%s: period=%.3f rotation=%.3f lattice=%.6f whitened=%.6f template=%.6f",
            path,
            score.selected_period,
            score.selected_rotation_degrees,
            score.joint_coherence,
            score.joint_whitened_match,
            score.canonical_template_score,
        )
    report_out.parent.mkdir(parents=True, exist_ok=True)
    report_out.write_text(
        json.dumps(
            {
                "schema_version": 13,
                "template": str(template_path),
                "periods": [float(value) for value in periods],
                "rotations_degrees": [float(value) for value in rotations],
                "input_angle": input_angle,
                "deskew_search": deskew_search,
                "input_scale": input_scale,
                "crop_fraction": crop_fraction,
                "jpeg_quality": jpeg_quality,
                "webp_quality": webp_quality,
                "input_rotation": int(input_rotation),
                "input_flip": input_flip,
                "orientation_search": orientation_search,
                "dihedral_search": dihedral_search,
                "registered_period": registered_period,
                "same_image_null": same_image_null,
                "patch_shift_consensus": patch_shift_consensus,
                "opponent_registered": opponent_registered,
                "amplitude_candidate_count": amplitude_candidate_count,
                "records": rows,
            },
            indent=2,
        )
        + "\n",
        encoding="utf-8",
    )
    log.info("Wrote %d lattice records: %s", len(rows), report_out)


if __name__ == "__main__":
    main()