Keep the pixel forensics in the library, drop the ai-score tooling

`scripts/ai_score.py` and the dataset scanner that fed it are gone: the detector
they trained is not something this project runs, and the corpus lived outside the
repository anyway. Nothing else referenced them.

The scanner's pixel layer was worth keeping, so it moves into the package as
`pixel_evidence.py` -- six families of scale-robust statistics (block-DCT histograms
and Benford deviation, FFT band energies and CFA peaks, high-pass residual, error
level, gradient, colour) measured in a single shared decode. The arithmetic was
verified against the scanner over 60 corpus images, families and artifacts alike,
before the scanner was removed; that comparison is no longer possible, which is why
the tests now pin behavior instead: determinism, empty-not-wrong on images too small
for a family, and one failing family not taking the others with it.

It has no consumer. Nothing in the package reads it, and the module says so.

`artifacts=True` returns the spatial layer -- perceptual hash, 128px thumbnail,
coarse ELA/residual/phase maps. Those identify the source image rather than describe
it, so they are opt-in and separate: everything else is a scalar or a fixed-length
histogram nothing can be reconstructed from.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

src/remove_ai_watermarks/pixel_evidence.py

@@ -0,0 +1,400 @@
+# pyright: reportUnknownMemberType=false, reportUnknownArgumentType=false, reportUnknownVariableType=false, reportMissingTypeStubs=false
+"""Experimental: the complete pixel-forensics layer for one image.
+
+STATUS
+
+Experimental and unused. Nothing in this package reads it -- not the provenance
+verdict, not removal, not the CLI. It is here because the research scanner that
+produced these measurements is gone, and the capability was worth keeping: whatever
+asks for pixel forensics next starts from a tested implementation instead of
+rebuilding one. Treat the shape as unstable until a caller exists.
+
+WHAT IS MEASURED
+
+One decode, then six families of scale-robust statistics over it:
+
+* ``dct`` -- AC coefficient histograms over the 8x8 block DCT, plus the deviation of
+  leading digits from Benford's law.
+* ``fft`` -- radial band energies of the log-magnitude spectrum, plus the
+  colour-filter-array periodicity peaks a demosaiced camera capture leaves.
+* ``noise`` -- standard deviation and kurtosis of a high-pass residual.
+* ``ela`` -- error level after a quality-90 JPEG re-save.
+* ``gradient`` -- gradient-magnitude histogram and Laplacian variance.
+* ``color`` -- 4x4x4 RGB histogram, mean saturation, mean value.
+
+and, in ``artifacts``, the spatial layer those statistics are computed from: a
+64-bit perceptual hash, a 128px JPEG thumbnail, and coarse ELA, noise-residual and
+FFT-phase maps.
+
+THE ARTIFACTS ARE NOT AGGREGATES
+
+Everything above ``artifacts`` is a scalar or a fixed-length histogram, and an image
+cannot be reconstructed from those. ``artifacts`` is different in kind: a thumbnail
+is a picture, a perceptual hash identifies one, and the coarse maps carry layout.
+Collecting them makes a record that identifies the source image, so a caller storing
+or forwarding them is handling image content, not statistics about it. That is why
+they are a separate field and not merged into the families.
+
+REQUIREMENTS
+
+Needs the ``pixels`` extra (numpy). Guard a call with :func:`is_available` when the
+caller must not hard-depend on it.
+"""
+
+from __future__ import annotations
+
+import base64
+import io
+import logging
+from dataclasses import dataclass, field
+from typing import TYPE_CHECKING, Any
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+logger = logging.getLogger(__name__)
+
+# Analysis resolution. Every statistic here is scale-robust, and a 2048px cap keeps
+# the FFT and the sliding-window residual bounded on a 100 MP input.
+MAX_SIDE = 2048
+# The eight lowest-frequency AC positions of the 8x8 block DCT, zig-zag order.
+AC_POSITIONS = ((0, 1), (1, 0), (1, 1), (0, 2), (2, 0), (2, 1), (1, 2), (0, 3))
+FFT_BANDS = 8
+# A Bayer CFA shows as symmetric peaks at half the Nyquist on the diagonals.
+BAYER_OFFSETS = ((1, 1), (1, -1))
+INSTALL_HINT = "install the pixel extra: uv add 'remove-ai-watermarks[pixels]'"
+
+
+@dataclass(frozen=True)
+class PixelEvidence:
+    """Pixel statistics for one image, and the spatial artifacts behind them.
+
+    ``decode`` carries the source dimensions, or ``{"error": ...}`` when the image
+    could not be decoded -- in which case every other field is empty. A family is also
+    empty when the image is too small for it (the block DCT needs 8x8, the FFT 32x32,
+    the residual 3x3), so a caller must treat every field as optional rather than
+    assume a fixed feature width.
+    """
+
+    path: Path
+    decode: dict[str, Any]
+    dct: dict[str, Any] = field(default_factory=dict[str, Any])
+    fft: dict[str, Any] = field(default_factory=dict[str, Any])
+    noise: dict[str, Any] = field(default_factory=dict[str, Any])
+    ela: dict[str, Any] = field(default_factory=dict[str, Any])
+    gradient: dict[str, Any] = field(default_factory=dict[str, Any])
+    color: dict[str, Any] = field(default_factory=dict[str, Any])
+    # Identifies the source image; see the module note. Empty unless asked for.
+    artifacts: dict[str, Any] = field(default_factory=dict[str, Any])
+
+    @property
+    def decoded(self) -> bool:
+        """False when the source could not be decoded at all."""
+        return "error" not in self.decode
+
+
+def is_available() -> bool:
+    """True when the optional pixel dependencies are installed."""
+    from remove_ai_watermarks.optional_deps import module_available
+
+    return module_available("numpy")
+
+
+def _numpy() -> Any:
+    from remove_ai_watermarks.optional_deps import module_available
+
+    if not module_available("numpy"):
+        raise RuntimeError(f"Pixel evidence needs numpy -- {INSTALL_HINT}")
+    import numpy as np
+
+    return np
+
+
+def _dct_matrix(np: Any, n: int = 8) -> Any:
+    """Orthonormal n x n DCT-II basis: M[i, j] = cos(pi (2j + 1) i / 2n)."""
+    i = np.arange(n)[:, None]
+    j = np.arange(n)[None, :]
+    m = np.cos(np.pi * (2 * j + 1) * i / (2 * n))
+    m[0, :] *= 1 / np.sqrt(2)
+    return m * np.sqrt(2 / n)
+
+
+def read_gray(image_path: Path) -> tuple[Any, Any, dict[str, Any]]:
+    """Decode to float32 grayscale (and RGB for colour stats), downscaled.
+
+    Pillow, not cv2, and the source dimensions are recorded BEFORE the downscale.
+    """
+    np = _numpy()
+    from PIL import Image
+
+    try:
+        with Image.open(image_path) as img:
+            info: dict[str, Any] = {"width": img.width, "height": img.height}
+            if max(img.size) > MAX_SIDE:
+                img.thumbnail((MAX_SIDE, MAX_SIDE), Image.Resampling.LANCZOS)
+            rgb = np.asarray(img.convert("RGB"), dtype=np.float32)
+            gray = np.asarray(img.convert("L"), dtype=np.float32)
+    except Exception as exc:
+        logger.debug("pixel decode failed for %s: %s", image_path, exc)
+        return None, None, {"error": f"{type(exc).__name__}: {exc}"}
+    return gray, rgb, info
+
+
+def dct_features(gray: Any) -> dict[str, Any]:
+    """AC coefficient histograms over the 8x8 block DCT + Benford deviation."""
+    np = _numpy()
+    height, width = gray.shape
+    h8, w8 = height // 8 * 8, width // 8 * 8
+    if h8 < 8 or w8 < 8:
+        return {}
+    basis = _dct_matrix(np)
+    bins = np.linspace(-20.5, 20.5, 22)
+    blocks = gray[:h8, :w8].reshape(h8 // 8, 8, w8 // 8, 8).swapaxes(1, 2)
+    coeff = np.einsum("ij,abjk,lk->abil", basis, blocks, basis)
+    hists = []
+    lead_vals: list[Any] = []
+    for dy, dx in AC_POSITIONS:
+        values = coeff[:, :, dy, dx].ravel()
+        hists.append(np.histogram(values, bins=bins)[0].tolist())
+        lead_vals.append(np.abs(values))
+    out: dict[str, Any] = {"dct_ac_hist": hists}
+    flat = np.abs(np.concatenate(lead_vals))
+    flat = flat[flat >= 1]
+    if flat.size > 100:
+        leading = (flat / 10 ** np.floor(np.log10(flat))).astype(int)
+        leading = leading[(leading >= 1) & (leading <= 9)]
+        if leading.size > 100:
+            observed = np.bincount(leading, minlength=10)[1:10] / leading.size
+            benford = np.log10(1 + 1 / np.arange(1, 10))
+            out["benford_mad"] = float(np.abs(observed - benford).mean())
+    return out
+
+
+def noise_residual_map(gray: Any) -> Any:
+    """High-pass residual, the map the noise statistics are computed from."""
+    np = _numpy()
+    from numpy.lib.stride_tricks import sliding_window_view
+
+    if gray.shape[0] < 3 or gray.shape[1] < 3:
+        return None
+    kernel = np.array([[-1.0, -1.0, -1.0], [-1.0, 8.0, -1.0], [-1.0, -1.0, -1.0]])
+    height, width = gray.shape
+    # kernel is float64, so the residual is float64 like the unchunked form
+    out = np.empty((height - 2, width - 2), dtype=np.float64)
+    # Row-chunked: the (window * kernel) temporary is ~150 MB at 2048px if
+    # materialized whole. Per-element 9-tap sums are computed in the same order,
+    # so the result is bit-identical to the unchunked form.
+    for y0 in range(0, height - 2, 256):
+        y1 = min(y0 + 256, height - 2)
+        window = sliding_window_view(gray[y0 : y1 + 2], (3, 3))
+        out[y0:y1] = (window * kernel).sum(axis=(-1, -2))
+    return out
+
+
+def noise_features(residual: Any) -> dict[str, Any]:
+    """High-pass residual std and kurtosis."""
+    flat = residual.ravel()
+    std = float(flat.std())
+    if std < 1e-9:
+        return {"noise_std": 0.0, "noise_kurtosis": 0.0}
+    z = (flat - flat.mean()) / std
+    return {"noise_std": std, "noise_kurtosis": float((z**4).mean() - 3.0)}
+
+
+def fft_decompose(gray: Any) -> tuple[Any, Any] | None:
+    """Log-magnitude (fftshifted) and phase of the image spectrum."""
+    np = _numpy()
+    if min(gray.shape) < 32:
+        return None
+    spectrum = np.fft.fftshift(np.fft.fft2(gray - gray.mean()))
+    return np.log1p(np.abs(spectrum)), np.angle(spectrum)
+
+
+def fft_features(mag: Any) -> dict[str, Any]:
+    """Radial magnitude band energies (no phase) + CFA periodicity peaks."""
+    np = _numpy()
+    height, width = mag.shape
+    cy, cx = height // 2, width // 2
+    # 1D broadcast instead of an mgrid: saves ~160 MB of int64 temporaries at
+    # 2048px. The squares are exact in float64 (values < 2^53), so band means
+    # are identical to the mgrid form.
+    r2y = (np.arange(height, dtype=np.float64) - cy) ** 2
+    r2x = (np.arange(width, dtype=np.float64) - cx) ** 2
+    radius = np.sqrt(r2y[:, None] + r2x[None, :])
+    r_max = radius.max()
+    bands = []
+    for index in range(FFT_BANDS):
+        mask = (radius >= r_max * index / FFT_BANDS) & (radius < r_max * (index + 1) / FFT_BANDS)
+        bands.append(float(mag[mask].mean()) if mask.any() else 0.0)
+    peaks = []
+    for dy, dx in BAYER_OFFSETS:
+        y, x = cy + dy * (height // 4), cx + dx * (width // 4)
+        neighborhood = mag[y - 2 : y + 3, x - 2 : x + 3]
+        peaks.append(float(neighborhood.max() - mag.mean()))
+    return {"fft_band_energy": bands, "cfa_peaks": peaks, "cfa_peak": max(peaks)}
+
+
+def ela_map(rgb: Any) -> Any:
+    """Absolute per-pixel error after a quality-90 JPEG re-save."""
+    np = _numpy()
+    from PIL import Image
+
+    try:
+        buffer = io.BytesIO()
+        Image.fromarray(rgb.astype(np.uint8)).save(buffer, "JPEG", quality=90)
+        buffer.seek(0)
+        resaved = np.asarray(Image.open(buffer).convert("RGB"), dtype=np.float32)
+    except Exception as exc:
+        logger.debug("ELA re-save failed: %s", exc)
+        return None
+    if resaved.shape != rgb.shape:
+        return None
+    return np.abs(rgb - resaved).mean(axis=-1)
+
+
+def ela_features(err: Any) -> dict[str, Any]:
+    """Error-level stats after a quality-90 JPEG re-save."""
+    np = _numpy()
+    return {"ela_mean": float(err.mean()), "ela_p95": float(np.percentile(err, 95))}
+
+
+def gradient_features(gray: Any) -> dict[str, Any]:
+    np = _numpy()
+    gy, gx = np.gradient(gray)
+    mag = np.sqrt(gx**2 + gy**2)
+    hist = np.histogram(mag, bins=10, range=(0, 255))[0].tolist()
+    laplacian = np.gradient(gy, axis=0) + np.gradient(gx, axis=1)
+    return {"gradient_hist": hist, "laplacian_var": float(laplacian.var())}
+
+
+def color_features(rgb: Any) -> dict[str, Any]:
+    np = _numpy()
+    small = rgb[::4, ::4]  # decimate; the histogram is position-blind anyway
+    bins = (small / 256 * 4).astype(int).clip(0, 3)
+    index = bins[..., 0] * 16 + bins[..., 1] * 4 + bins[..., 2]
+    hist = np.bincount(index.ravel(), minlength=64).tolist()
+    mx = small.max(axis=-1)
+    mn = small.min(axis=-1)
+    saturation = np.where(mx > 0, (mx - mn) / np.maximum(mx, 1e-6), 0)
+    return {
+        "color_hist_4x4x4": hist,
+        "saturation_mean": float(saturation.mean()),
+        "value_mean": float(mx.mean() / 255),
+    }
+
+
+def perceptual_hash(gray: Any) -> str:
+    """64-bit DCT perceptual hash. Identifies an image; see the module note."""
+    np = _numpy()
+    from PIL import Image
+
+    small = np.asarray(Image.fromarray(gray.astype(np.float32), mode="F").resize((32, 32), Image.Resampling.LANCZOS))
+    basis = _dct_matrix(np, 32)
+    coeff = basis @ small @ basis.T
+    low = coeff[:8, :8].ravel()[1:]  # drop DC
+    bits = low > np.median(low)
+    return f"{int(''.join('1' if bit else '0' for bit in bits), 2):016x}"
+
+
+def _coarse(np: Any, arr: Any, side: int = 64) -> Any:
+    """Downscale a 2D map to at most ``side`` on the long edge."""
+    from PIL import Image
+
+    height, width = arr.shape
+    if max(height, width) <= side:
+        return arr
+    img = Image.fromarray(arr.astype(np.float32), mode="F")
+    img.thumbnail((side, side), Image.Resampling.BILINEAR)
+    return np.asarray(img)
+
+
+def _array_payload(arr: Any) -> dict[str, Any]:
+    return {
+        "shape": list(arr.shape),
+        "dtype": str(arr.dtype),
+        "base64": base64.b64encode(arr.tobytes()).decode("ascii"),
+    }
+
+
+def spatial_artifacts(gray: Any, rgb: Any, *, ela: Any, residual: Any, phase: Any) -> dict[str, Any]:
+    """Perceptual hash, thumbnail, and coarse ELA / residual / phase maps.
+
+    These identify the source image rather than describe it -- see the module note.
+    The maps are the ones the statistics were computed from, passed in rather than
+    recomputed.
+    """
+    np = _numpy()
+    from PIL import Image
+
+    out: dict[str, Any] = {"phash": perceptual_hash(gray)}
+    thumbnail = Image.fromarray(rgb.astype(np.uint8))
+    thumbnail.thumbnail((128, 128), Image.Resampling.LANCZOS)
+    buffer = io.BytesIO()
+    thumbnail.save(buffer, "JPEG", quality=70)
+    out["thumbnail_jpeg_b64"] = base64.b64encode(buffer.getvalue()).decode("ascii")
+
+    if ela is not None:
+        out["ela_map"] = _array_payload(_coarse(np, ela))
+    if residual is not None:
+        clipped = np.clip(residual / 4.0, -1, 1)
+        out["noise_residual"] = _array_payload(_coarse(np, (clipped * 127).astype(np.int8)))
+    if phase is not None:
+        out["fft_phase"] = _array_payload(_coarse(np, phase.astype(np.float32), 32))
+    return out
+
+
+def extract_pixel_evidence(image_path: Path, *, artifacts: bool = False) -> PixelEvidence:
+    """Measure every pixel-statistic family for one image in a single decode.
+
+    The image is decoded ONCE and the intermediate maps (high-pass residual, ELA
+    error, FFT magnitude and phase) are computed once and shared, because the
+    residual's sliding window and the ELA re-save are the two expensive steps and
+    each family would otherwise redo them.
+
+    A family that fails or does not apply is left empty rather than raising: an
+    undecodable file, or one too small for the block DCT, still returns a
+    :class:`PixelEvidence` whose ``decoded`` / empty fields say so. Missing numpy is
+    the one hard error, since then nothing can be measured at all.
+
+    Args:
+        image_path: Path to the image. Any container Pillow can open.
+        artifacts: Also return the spatial layer -- perceptual hash, thumbnail and
+            coarse maps. Off by default: those identify the source image, so asking
+            for them is a decision the caller makes explicitly.
+
+    Returns:
+        A :class:`PixelEvidence`.
+    """
+    gray, rgb, info = read_gray(image_path)
+    if gray is None or rgb is None:
+        return PixelEvidence(path=image_path, decode=info)
+
+    residual = noise_residual_map(gray)
+    spectrum = fft_decompose(gray)
+    error = ela_map(rgb)
+
+    families: dict[str, dict[str, Any]] = {}
+    for name, compute in (
+        ("noise", lambda: noise_features(residual) if residual is not None else {}),
+        ("fft", lambda: fft_features(spectrum[0]) if spectrum is not None else {}),
+        ("ela", lambda: ela_features(error) if error is not None else {}),
+        ("dct", lambda: dct_features(gray)),
+        ("gradient", lambda: gradient_features(gray)),
+        ("color", lambda: color_features(rgb)),
+    ):
+        try:
+            families[name] = compute()
+        except Exception as exc:  # one bad family must not lose the other five
+            logger.debug("pixel family %s failed for %s: %s", name, image_path, exc)
+            families[name] = {"error": f"{type(exc).__name__}: {exc}"}
+
+    if artifacts:
+        try:
+            families["artifacts"] = spatial_artifacts(
+                gray, rgb, ela=error, residual=residual, phase=spectrum[1] if spectrum is not None else None
+            )
+        except Exception as exc:
+            logger.debug("pixel artifacts failed for %s: %s", image_path, exc)
+            families["artifacts"] = {"error": f"{type(exc).__name__}: {exc}"}
+
+    return PixelEvidence(path=image_path, decode=info, **families)