# rtt_jwst_sidecar.py
from __future__ import annotations
 
import json
import argparse
from dataclasses import dataclass
from pathlib import Path
from typing import Dict, Any, Optional
 
import numpy as np
from jwst import datamodels
 
 
@dataclass
class Thresholds:
    # Keep these small and boring. Expand later per instrument.
    max_nan_frac: float = 0.01
    max_background_sigma: float = 5.0
    min_snr_proxy: float = 2.0
 
 
def robust_stats(x: np.ndarray) -> Dict[str, float]:
    x = x[np.isfinite(x)]
    if x.size == 0:
        return {"median": float("nan"), "mad": float("nan"), "sigma": float("nan")}
    med = float(np.median(x))
    mad = float(np.median(np.abs(x - med)))
    sigma = float(1.4826 * mad)  # MAD->sigma proxy
    return {"median": med, "mad": mad, "sigma": sigma}
 
 
def compute_rtt_signals(model) -> Dict[str, float]:
    """
    RTT signals (bounded, minimal):
    - noise: robust sigma of background-ish pixels
    - stability: inverse of NaN fraction + noise inflation
    - internal: consistency of error array vs science array (proxy)
    - memory: not true 'memory' yet; use header completeness proxy
    """
    data = getattr(model, "data", None)
    err = getattr(model, "err", None)
 
    if data is None:
        return {"noise": 1.0, "stability": 0.0, "internal": 0.0, "memory": 0.0}
 
    data = np.array(data, dtype=float)
    nan_frac = float(np.mean(~np.isfinite(data)))
 
    # crude background proxy: use lower half of pixel distribution
    finite = data[np.isfinite(data)]
    if finite.size == 0:
        bg_stats = {"sigma": float("nan")}
    else:
        cutoff = np.percentile(finite, 50)
        bg = finite[finite <= cutoff]
        bg_stats = robust_stats(bg)
 
    noise = float(np.clip(bg_stats["sigma"] / (np.nanstd(finite) + 1e-12), 0.0, 1.0)) if finite.size else 1.0
 
    # internal consistency proxy: err should not be missing and should correlate with |data| a bit
    internal = 0.0
    if err is not None:
        err = np.array(err, dtype=float)
        ok = np.isfinite(err) & np.isfinite(data)
        if np.any(ok):
            a = np.abs(data[ok]).ravel()
            b = np.abs(err[ok]).ravel()
            if a.size > 50:
                corr = np.corrcoef(a[:5000], b[:5000])[0, 1]
                internal = float(np.clip((corr + 1) / 2, 0.0, 1.0))
 
    # memory proxy: do we have key metadata for downstream traceability?
    meta = getattr(model, "meta", None)
    memory = 0.0
    if meta is not None:
        keys = [
            ("observation", "program_number"),
            ("observation", "observation_number"),
            ("instrument", "name"),
            ("exposure", "type"),
        ]
        present = 0
        for a, b in keys:
            v = getattr(getattr(meta, a, None), b, None)
            present += int(v is not None)
        memory = present / len(keys)
 
    # stability: penalize NaNs + high noise
    stability = float(np.clip(1.0 - (0.7 * nan_frac + 0.3 * noise), 0.0, 1.0))
 
    return {
        "nan_frac": nan_frac,
        "bg_sigma_proxy": float(bg_stats["sigma"]),
        "noise": noise,
        "internal": internal,
        "memory": memory,
        "stability": stability,
    }
 
 
def classify_zone(sig: Dict[str, float], thr: Thresholds) -> str:
    # Simple and explicit. Expand later.
    if sig["nan_frac"] > thr.max_nan_frac:
        return "Transit Verge"
    if sig["bg_sigma_proxy"] > thr.max_background_sigma:
        return "Transit Verge"
    if sig["stability"] > 0.8 and sig["internal"] > 0.6:
        return "Lagrange Calm"
    if sig["noise"] < 0.3 and sig["memory"] > 0.7:
        return "Deep Quiet"
    return "Echo Belt"
 
 
def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("jwst_fits", help="Input JWST product (rate/cal/i2d/s3d/x1d, etc.)")
    ap.add_argument("--out", default=None, help="Output JSON path (default: <input>.rtt.json)")
    args = ap.parse_args()
 
    in_path = Path(args.jwst_fits)
    out_path = Path(args.out) if args.out else in_path.with_suffix(in_path.suffix + ".rtt.json")
 
    thr = Thresholds()
    with datamodels.open(str(in_path)) as model:
        sig = compute_rtt_signals(model)
        zone = classify_zone(sig, thr)
 
        report: Dict[str, Any] = {
            "rtt_version": "0.1",
            "input": str(in_path),
            "zone": zone,
            "signals": sig,
            "thresholds": thr.__dict__,
            "notes": "Sidecar QA only; does not modify calibrated arrays.",
        }
 
    out_path.write_text(json.dumps(report, indent=2))
    print(f"Wrote {out_path}")
 
 
if __name__ == "__main__":
    main()

from dataclasses import dataclass
 
@dataclass(frozen=True)
class NircamImagingThresholds:
    # data integrity
    max_nan_frac: float = 0.01
    max_do_not_use_frac: float = 0.02
 
    # background / noise
    max_bg_sigma_e_per_s: float = 8.0          # conservative starter; tune per filter
    max_bg_gradient_e_per_s_per_pix: float = 0.15
 
    # sharpness proxy (dimensionless)
    min_sharpness: float = 0.010               # too low => smear/over-smoothing
    max_sharpness: float = 0.120               # too high => unmasked CRs/artifacts
 
    # zone thresholds
    calm_stability_min: float = 0.80
    verge_risk_min: float = 0.65

from __future__ import annotations
 
import argparse
import json
from dataclasses import asdict, dataclass
from pathlib import Path
from typing import Any, Dict, Optional, Tuple
 
import numpy as np
from jwst import datamodels
 
from profiles_nircam_imaging import NircamImagingThresholds
 
 
def clamp(x: float, lo: float = 0.0, hi: float = 1.0) -> float:
    return max(lo, min(hi, x))
 
 
def robust_sigma(x: np.ndarray) -> float:
    x = x[np.isfinite(x)]
    if x.size == 0:
        return float("nan")
    med = np.median(x)
    mad = np.median(np.abs(x - med))
    return float(1.4826 * mad)
 
 
def estimate_background(data: np.ndarray) -> Tuple[float, float]:
    """
    Returns: (bg_level, bg_sigma) using a conservative lower-half estimate.
    """
    finite = data[np.isfinite(data)]
    if finite.size == 0:
        return float("nan"), float("nan")
    cutoff = np.percentile(finite, 50)
    bg = finite[finite <= cutoff]
    return float(np.median(bg)), robust_sigma(bg)
 
 
def gradient_proxy(data: np.ndarray) -> float:
    """
    Simple gradient magnitude proxy: median(|dx| + |dy|) on a decimated grid.
    """
    x = data[::4, ::4]
    x = np.where(np.isfinite(x), x, np.nan)
    dx = np.nanmedian(np.abs(np.diff(x, axis=1)))
    dy = np.nanmedian(np.abs(np.diff(x, axis=0)))
    if not np.isfinite(dx) and not np.isfinite(dy):
        return float("nan")
    return float(np.nan_to_num(dx, nan=0.0) + np.nan_to_num(dy, nan=0.0))
 
 
def sharpness_proxy(data: np.ndarray) -> float:
    """
    High-frequency energy proxy:
    - compute simple Laplacian magnitude on decimated data
    - normalize by robust sigma to make it somewhat scale-invariant
    """
    x = data[::2, ::2]
    x = np.where(np.isfinite(x), x, 0.0)
 
    lap = (
        -4.0 * x
        + np.roll(x, 1, 0) + np.roll(x, -1, 0)
        + np.roll(x, 1, 1) + np.roll(x, -1, 1)
    )
    hf = np.median(np.abs(lap))
    sig = robust_sigma(x)
    if not np.isfinite(sig) or sig <= 0:
        return float("nan")
    return float(hf / sig)
 
 
def dq_fractions(dq: Optional[np.ndarray]) -> Dict[str, float]:
    """
    Without importing jwst.datamodels.dqflags, we keep it minimal:
    - fraction nonzero (any flag) as a coarse integrity signal
    If you want specific bits (DO_NOT_USE/SATURATED), we can add dqflags later.
    """
    if dq is None:
        return {"dq_nonzero_frac": float("nan")}
    dq = np.array(dq)
    total = dq.size
    if total == 0:
        return {"dq_nonzero_frac": float("nan")}
    return {"dq_nonzero_frac": float(np.count_nonzero(dq) / total)}
 
 
def compute_signals(model) -> Dict[str, float]:
    data = getattr(model, "data", None)
    dq = getattr(model, "dq", None)
 
    if data is None:
        return {
            "nan_frac": 1.0,
            "dq_nonzero_frac": float("nan"),
            "bg_sigma": float("nan"),
            "bg_gradient": float("nan"),
            "sharpness": float("nan"),
            "wcs_present": 0.0,
        }
 
    data = np.array(data, dtype=float)
    nan_frac = float(np.mean(~np.isfinite(data)))
 
    bg_level, bg_sigma = estimate_background(data)
    bg_grad = gradient_proxy(data)
    sharp = sharpness_proxy(data)
 
    dq_info = dq_fractions(dq)
 
    # WCS presence: Stage 2+ should have it in metadata products
    wcs_present = 0.0
    try:
        wcs_present = 1.0 if getattr(model.meta, "wcs", None) is not None else 0.0
    except Exception:
        wcs_present = 0.0
 
    return {
        "nan_frac": nan_frac,
        **dq_info,
        "bg_level": float(bg_level),
        "bg_sigma": float(bg_sigma),
        "bg_gradient": float(bg_grad),
        "sharpness": float(sharp),
        "wcs_present": float(wcs_present),
    }
 
 
def rtt_scores(sig: Dict[str, float], thr: NircamImagingThresholds) -> Dict[str, float]:
    # Noise: map bg_sigma into 0..1 where 1 is “too noisy”
    if np.isfinite(sig["bg_sigma"]):
        noise = clamp(sig["bg_sigma"] / max(1e-6, thr.max_bg_sigma_e_per_s))
    else:
        noise = 1.0
 
    # Risk: integrity + WCS missing + sharpness out-of-family
    dq_nonzero = sig.get("dq_nonzero_frac", float("nan"))
    dq_term = 0.0 if not np.isfinite(dq_nonzero) else clamp(dq_nonzero / max(1e-6, thr.max_do_not_use_frac))
    nan_term = clamp(sig["nan_frac"] / max(1e-6, thr.max_nan_frac))
    wcs_term = 1.0 - sig.get("wcs_present", 0.0)
 
    sharp = sig.get("sharpness", float("nan"))
    sharp_term = 0.0
    if np.isfinite(sharp):
        if sharp < thr.min_sharpness:
            sharp_term = clamp((thr.min_sharpness - sharp) / thr.min_sharpness)
        elif sharp > thr.max_sharpness:
            sharp_term = clamp((sharp - thr.max_sharpness) / thr.max_sharpness)
 
    risk = clamp(0.35 * dq_term + 0.35 * nan_term + 0.20 * wcs_term + 0.10 * sharp_term)
 
    # Stability: high when noise/risk/gradient are low
    grad = sig.get("bg_gradient", float("nan"))
    grad_term = 0.0 if not np.isfinite(grad) else clamp(grad / max(1e-6, thr.max_bg_gradient_e_per_s_per_pix))
    stability = clamp(1.0 - (0.45 * noise + 0.40 * risk + 0.15 * grad_term))
 
    # Internal: consistency proxy (here: WCS present + finite background stats)
    internal = clamp(0.5 * sig.get("wcs_present", 0.0) + 0.5 * float(np.isfinite(sig.get("bg_sigma", float("nan")))))
 
    # Memory: we keep “memory” as metadata completeness proxy
    memory = internal
 
    coherence = clamp(0.40 * stability + 0.25 * internal + 0.15 * memory - 0.10 * noise - 0.10 * risk)
 
    return {
        "noise": float(noise),
        "risk": float(risk),
        "stability": float(stability),
        "internal": float(internal),
        "memory": float(memory),
        "coherence": float(coherence),
    }
 
 
def classify_zone(scores: Dict[str, float], thr: NircamImagingThresholds) -> Tuple[str, str]:
    # Explain “why” in plain terms
    if scores["risk"] >= thr.verge_risk_min:
        return "Transit Verge", "Risk high: integrity/WCS/sharpness boundary crossed"
    if scores["stability"] >= thr.calm_stability_min and scores["risk"] < 0.30:
        return "Lagrange Calm", "Stable + low risk"
    if scores["noise"] < 0.30 and scores["internal"] > 0.70:
        return "Deep Quiet", "Clean background + metadata coherence"
    return "Echo Belt", "Nominal but not fully calm (trend watch)"
 
 
def main() -> None:
    ap = argparse.ArgumentParser(description="RTT sidecar QA for JWST NIRCam imaging products.")
    ap.add_argument("fits", help="Input JWST product (rate/cal/i2d).")
    ap.add_argument("--out", default=None, help="Output JSON path (default: <input>.rtt.json).")
    ap.add_argument("--thr", default=None, help="Optional thresholds JSON (override defaults).")
    args = ap.parse_args()
 
    in_path = Path(args.fits)
    out_path = Path(args.out) if args.out else in_path.with_suffix(in_path.suffix + ".rtt.json")
 
    thr = NircamImagingThresholds()
    if args.thr:
        override = json.loads(Path(args.thr).read_text(encoding="utf-8"))
        thr = NircamImagingThresholds(**{**asdict(thr), **override})
 
    with datamodels.open(str(in_path)) as model:
        sig = compute_signals(model)
        scores = rtt_scores(sig, thr)
        zone, why = classify_zone(scores, thr)
 
        meta = {}
        try:
            meta = {
                "instrument": getattr(model.meta.instrument, "name", None),
                "exposure_type": getattr(model.meta.exposure, "type", None),
                "detector": getattr(model.meta.instrument, "detector", None),
                "filter": getattr(model.meta.instrument, "filter", None),
                "pupil": getattr(model.meta.instrument, "pupil", None),
                "pipeline_stage_hint": in_path.name.split("_")[-1],
            }
        except Exception:
            pass
 
    report: Dict[str, Any] = {
        "rtt_version": "0.1-nircam-imaging",
        "input": str(in_path),
        "zone": zone,
        "why": why,
        "signals": sig,
        "scores": scores,
        "thresholds": asdict(thr),
        "meta": meta,
        "notes": "Sidecar QA only; does not modify calibrated arrays. Designed to run post Stage 2/3 products.",
    }
 
    out_path.write_text(json.dumps(report, indent=2), encoding="utf-8")
    print(f"Wrote {out_path}")
 
 
if __name__ == "__main__":
    main()

from dataclasses import dataclass
 
@dataclass(frozen=True)
class NircamMosaicThresholds:
    # --- Seam coherence ---
    max_seam_contrast: float              # seam residual / local bg sigma
    max_seam_abs_e_per_s: float           # absolute seam jump
 
    # --- Background uniformity ---
    max_tile_bg_sigma_ratio: float        # tile bg sigma / global bg sigma
    max_global_gradient_e_per_s_per_pix: float
 
    # --- Sharpness continuity ---
    min_sharpness: float
    max_sharpness: float
    max_sharpness_jump: float             # across seam
 
    # --- Astrometric plausibility ---
    max_wcs_offset_arcsec: float
    max_pixel_scale_frac: float
 
    # --- Zone thresholds ---
    calm_stability_min: float
    verge_risk_min: float
 
 
# -------------------------
# NIRCam Short-Wavelength
# -------------------------
 
NIRCAM_SW_MOSAIC = NircamMosaicThresholds(
    # seams
    max_seam_contrast=1.5,
    max_seam_abs_e_per_s=0.8,
 
    # background
    max_tile_bg_sigma_ratio=1.25,
    max_global_gradient_e_per_s_per_pix=0.10,
 
    # sharpness
    min_sharpness=0.012,
    max_sharpness=0.120,
    max_sharpness_jump=0.020,
 
    # astrometry
    max_wcs_offset_arcsec=0.05,
    max_pixel_scale_frac=0.01,
 
    # zones
    calm_stability_min=0.82,
    verge_risk_min=0.65,
)
 
 
# -------------------------
# NIRCam Long-Wavelength
# -------------------------
 
NIRCAM_LW_MOSAIC = NircamMosaicThresholds(
    # seams (LW tolerates slightly larger residuals)
    max_seam_contrast=1.8,
    max_seam_abs_e_per_s=1.2,
 
    # background
    max_tile_bg_sigma_ratio=1.35,
    max_global_gradient_e_per_s_per_pix=0.15,
 
    # sharpness
    min_sharpness=0.008,
    max_sharpness=0.090,
    max_sharpness_jump=0.025,
 
    # astrometry
    max_wcs_offset_arcsec=0.07,
    max_pixel_scale_frac=0.015,
 
    # zones
    calm_stability_min=0.78,
    verge_risk_min=0.65,
)

#!/usr/bin/env python3
"""
rtt_jwst_nircam_i2d_sidecar.py
==============================
 
RTT sidecar QA for JWST NIRCam mosaic-level products (Stage 3), focused on:
- seam/coadd coherence for public-facing trust
- background uniformity / gradients
- lightweight WCS plausibility checks
- optional seam overlay PNGs (OFF by default)
 
This tool is READ-ONLY on input FITS. It emits:
- <input>.rtt.json (full report)
- <input>.rtt.csv  (single-row summary for batch dashboards)
- optional <input>.rtt_seams.png (visual overlay for quick human verification)
 
Intended inputs:
- NIRCam Stage 3 mosaic: *i2d.fits (ImageModel-like)
Also works on other 2D combined images if they follow JWST datamodel conventions.
 
Dependencies:
- jwst (stpipe/datamodels)
- numpy
- (optional) matplotlib for PNG overlays
 
Usage:
  python rtt_jwst_nircam_i2d_sidecar.py path/to/file_i2d.fits --band SW
  python rtt_jwst_nircam_i2d_sidecar.py path/to/file_i2d.fits --band LW --png
 
Batch usage:
  python rtt_jwst_nircam_i2d_sidecar.py *.fits --band SW --csv-out logs/rtt_i2d_summary.csv
 
Notes:
- "Tile segmentation" here is derived from drizzle weight maps when available.
  If weights are absent, we fall back to a coarse grid segmentation.
- Seam detection is computed on the weight-boundary edges and evaluated in a local band.
- WCS checks are lightweight and intentionally non-invasive (no astrometric solving).
"""
 
from __future__ import annotations
 
import argparse
import csv
import json
import math
from dataclasses import asdict
from pathlib import Path
from typing import Any, Dict, Optional, Tuple
 
import numpy as np
from jwst import datamodels
 
from profiles_nircam_sw_lw import NIRCAM_SW_MOSAIC, NIRCAM_LW_MOSAIC, NircamMosaicThresholds
 
 
# ----------------------------
# Utilities
# ----------------------------
 
def clamp(x: float, lo: float = 0.0, hi: float = 1.0) -> float:
    return max(lo, min(hi, x))
 
 
def finite(x: np.ndarray) -> np.ndarray:
    return np.isfinite(x)
 
 
def robust_median(x: np.ndarray) -> float:
    x = x[np.isfinite(x)]
    if x.size == 0:
        return float("nan")
    return float(np.median(x))
 
 
def robust_sigma(x: np.ndarray) -> float:
    """
    Robust sigma via MAD, stable under outliers.
    """
    x = x[np.isfinite(x)]
    if x.size == 0:
        return float("nan")
    med = np.median(x)
    mad = np.median(np.abs(x - med))
    return float(1.4826 * mad)
 
 
def downsample2d(a: np.ndarray, step: int) -> np.ndarray:
    if step <= 1:
        return a
    return a[::step, ::step]
 
 
def safe_float(v: Any) -> Optional[float]:
    try:
        if v is None:
            return None
        return float(v)
    except Exception:
        return None
 
 
# ----------------------------
# Input interpretation
# ----------------------------
 
def pick_thresholds(band: str) -> NircamMosaicThresholds:
    b = band.strip().upper()
    if b == "SW":
        return NIRCAM_SW_MOSAIC
    if b == "LW":
        return NIRCAM_LW_MOSAIC
    raise ValueError("band must be SW or LW")
 
 
def get_arrays(model) -> Tuple[np.ndarray, Optional[np.ndarray]]:
    """
    Returns (data, weight) arrays.
    Weight might be absent depending on model/product.
    """
    data = getattr(model, "data", None)
    if data is None:
        raise ValueError("Input model has no 'data' array")
 
    data = np.array(data, dtype=float)
 
    w = getattr(model, "wht", None)
    if w is None:
        w = getattr(model, "weight", None)  # some combined products
    if w is not None:
        w = np.array(w, dtype=float)
 
    return data, w
 
 
def infer_band_from_meta(model) -> Optional[str]:
    """
    Try to infer SW/LW band from filter/wavelength hints.
    Conservative: return None if unsure.
    """
    try:
        filt = getattr(model.meta.instrument, "filter", None)
        if filt and isinstance(filt, str):
            # Common NIRCam filters: F070W...F200W (SW), F277W...F444W (LW)
            # This is heuristic; users can override with --band.
            if filt.startswith("F") and filt.endswith(("W", "M", "N")) and len(filt) >= 5:
                n = int(filt[1:4])
                if n <= 210:
                    return "SW"
                if n >= 250:
                    return "LW"
    except Exception:
        pass
    return None
 
 
# ----------------------------
# Tile segmentation from weights
# ----------------------------
 
def segmentation_from_weights(wht: Optional[np.ndarray], grid: int = 8) -> Tuple[np.ndarray, Dict[str, Any]]:
    """
    Produce a coarse "tile id" map (same shape as data) and metadata describing it.
 
    Preferred: derive tile edges from weight discontinuities (drizzle footprint).
    Fallback: uniform grid segmentation.
    """
    info: Dict[str, Any] = {"method": None, "grid": grid}
 
    if wht is None or wht.size == 0 or not np.any(np.isfinite(wht)):
        info["method"] = "grid_fallback_no_weight"
        return grid_segmentation((0, 0), grid, shape=None, wht=None)
 
    # Normalize and downsample to find structure without overcost
    h, w = wht.shape
    step = 4 if max(h, w) > 2000 else 2
    ws = downsample2d(np.where(np.isfinite(wht), wht, 0.0), step)
 
    # Identify "coverage" region
    cov = ws > 0
    if np.count_nonzero(cov) < 1000:
        info["method"] = "grid_fallback_sparse_weight"
        return grid_segmentation((h, w), grid, shape=(h, w), wht=wht)
 
    # Weight gradient magnitude as proxy for boundaries
    gx = np.abs(np.diff(ws, axis=1, prepend=ws[:, :1]))
    gy = np.abs(np.diff(ws, axis=0, prepend=ws[:1, :]))
    g = gx + gy
 
    # Threshold boundaries relative to robust scale
    g_med = np.median(g[cov]) if np.any(cov) else 0.0
    g_sig = robust_sigma(g[cov]) if np.any(cov) else 0.0
    thr = g_med + 4.0 * (g_sig if np.isfinite(g_sig) else 0.0)
 
    boundary = (g > thr) & cov
    boundary_frac = float(np.count_nonzero(boundary) / max(1, np.count_nonzero(cov)))
    info["method"] = "weight_boundary"
    info["downsample_step"] = step
    info["boundary_threshold"] = float(thr)
    info["boundary_frac"] = boundary_frac
 
    # Convert boundary map to full resolution seam mask by nearest expansion
    boundary_full = np.zeros((h, w), dtype=bool)
    boundary_full[::step, ::step] = boundary
    # Fill gaps by simple dilation-like expansion (cheap)
    boundary_full = expand_mask(boundary_full, radius=step)
 
    # Tile segmentation itself: for v0.1 mosaic we use grid segments, but we keep
    # boundary_full for seam evaluation. This avoids heavy connected-component tiling.
    tile_id = grid_segmentation((h, w), grid, shape=(h, w), wht=wht)[0]
    info["tile_segmentation"] = "grid"
    return tile_id, {"tile": info, "seam_mask_method": "weight_boundary+grid_tiles", "boundary_frac": boundary_frac}
 
 
def grid_segmentation(size: Tuple[int, int], grid: int, shape: Optional[Tuple[int, int]], wht: Optional[np.ndarray]):
    """
    Create uniform grid tile labels over the valid footprint.
    """
    if shape is None:
        # caller passed nonsense; just return placeholder
        return np.zeros((1, 1), dtype=int), {"method": "grid_placeholder"}
 
    h, w = shape
    tile = np.full((h, w), -1, dtype=int)
 
    # Valid footprint based on weight>0 if present else finite data
    if wht is not None:
        valid = np.isfinite(wht) & (wht > 0)
    else:
        valid = np.ones((h, w), dtype=bool)
 
    gh = max(1, grid)
    gw = max(1, grid)
    dh = h / gh
    dw = w / gw
 
    tid = 0
    for i in range(gh):
        y0 = int(i * dh)
        y1 = int((i + 1) * dh) if i < gh - 1 else h
        for j in range(gw):
            x0 = int(j * dw)
            x1 = int((j + 1) * dw) if j < gw - 1 else w
            m = valid[y0:y1, x0:x1]
            if np.count_nonzero(m) > 0:
                tile[y0:y1, x0:x1][m] = tid
            tid += 1
 
    return tile, {"method": "grid", "grid": grid, "tiles": gh * gw}
 
 
def expand_mask(m: np.ndarray, radius: int = 2) -> np.ndarray:
    """
    Cheap expansion (approx dilation) without scipy.
    """
    r = max(1, int(radius))
    out = m.copy()
    for dy in range(-r, r + 1):
        for dx in range(-r, r + 1):
            if dy == 0 and dx == 0:
                continue
            out |= np.roll(np.roll(m, dy, axis=0), dx, axis=1)
    return out
 
 
# ----------------------------
# Seam detection & metrics
# ----------------------------
 
def seam_mask_from_tile_ids(tile_id: np.ndarray, valid: np.ndarray) -> np.ndarray:
    """
    Seam mask where adjacent pixels belong to different tile IDs.
    """
    seam = np.zeros_like(valid, dtype=bool)
    t = tile_id
 
    # Horizontal seams
    seam[:, 1:] |= (t[:, 1:] != t[:, :-1]) & valid[:, 1:] & valid[:, :-1]
    # Vertical seams
    seam[1:, :] |= (t[1:, :] != t[:-1, :]) & valid[1:, :] & valid[:-1, :]
 
    return seam
 
 
def seam_metrics(
    data: np.ndarray,
    wht: Optional[np.ndarray],
    tile_id: np.ndarray,
    seam_mask: np.ndarray,
    band_width: int = 8,
) -> Dict[str, float]:
    """
    Compute seam coherence metrics:
    - seam_abs: median(|local_mean_left - local_mean_right|) approximated via seam-band residuals
    - seam_contrast: seam_abs / local_bg_sigma
    - seam_density: fraction of valid pixels near seams
    """
    if wht is not None:
        valid = np.isfinite(data) & np.isfinite(wht) & (wht > 0)
    else:
        valid = np.isfinite(data)
 
    if np.count_nonzero(valid) < 1000:
        return {
            "seam_abs_e_per_s": float("nan"),
            "seam_contrast": float("nan"),
            "seam_density": float("nan"),
        }
 
    # Build a seam band mask
    band = expand_mask(seam_mask & valid, radius=band_width)
    band &= valid
 
    seam_density = float(np.count_nonzero(band) / np.count_nonzero(valid))
 
    if np.count_nonzero(band) < 500:
        return {
            "seam_abs_e_per_s": float("nan"),
            "seam_contrast": float("nan"),
            "seam_density": seam_density,
        }
 
    # Residual relative to local background (use robust median from non-bright pixels)
    # Use a conservative background sample: lower half of pixel distribution within valid.
    v = data[valid]
    cutoff = np.percentile(v, 50)
    bg = v[v <= cutoff]
    bg_sigma = robust_sigma(bg)
 
    # seam_abs: robust sigma of differences between band pixels and global median
    # We treat "seaminess" as local deviation structure near seams.
    bg_med = float(np.median(bg)) if bg.size else float(np.median(v))
    seam_dev = data[band] - bg_med
    seam_abs = float(np.median(np.abs(seam_dev)))
 
    seam_contrast = float(seam_abs / bg_sigma) if np.isfinite(bg_sigma) and bg_sigma > 0 else float("nan")
 
    return {
        "seam_abs_e_per_s": seam_abs,
        "seam_contrast": seam_contrast,
        "seam_density": seam_density,
        "bg_sigma_e_per_s": float(bg_sigma),
        "bg_med_e_per_s": float(bg_med),
    }
 
 
# ----------------------------
# Mosaic uniformity metrics
# ----------------------------
 
def tile_stats(data: np.ndarray, wht: Optional[np.ndarray], tile_id: np.ndarray) -> Dict[str, Any]:
    """
    Compute tile-to-tile background sigma ratios and background levels on a coarse tile segmentation.
    """
    if wht is not None:
        valid = np.isfinite(data) & np.isfinite(wht) & (wht > 0)
    else:
        valid = np.isfinite(data)
 
    v = data[valid]
    if v.size < 1000:
        return {"tile_bg_sigma_ratio": float("nan"), "tile_bg_med_spread": float("nan")}
 
    # Global background sample = lower half
    cutoff = np.percentile(v, 50)
    bg = v[v <= cutoff]
    g_sig = robust_sigma(bg)
    g_med = float(np.median(bg))
 
    # Per tile background sigma, use a minimum sample requirement
    tile_ids = np.unique(tile_id[valid])
    tile_ids = tile_ids[tile_ids >= 0]
 
    sigs = []
    meds = []
    for tid in tile_ids:
        m = valid & (tile_id == tid)
        x = data[m]
        if x.size < 2000:
            continue
        c = np.percentile(x, 50)
        xb = x[x <= c]
        if xb.size < 500:
            continue
        sigs.append(robust_sigma(xb))
        meds.append(float(np.median(xb)))
 
    if len(sigs) < 3 or not np.isfinite(g_sig) or g_sig <= 0:
        return {
            "tile_bg_sigma_ratio": float("nan"),
            "tile_bg_med_spread": float("nan"),
            "global_bg_sigma_e_per_s": float(g_sig),
            "global_bg_med_e_per_s": float(g_med),
            "tiles_used": int(len(sigs)),
        }
 
    sigs = np.array(sigs, dtype=float)
    meds = np.array(meds, dtype=float)
 
    ratio = float(np.nanmax(sigs) / np.nanmedian(sigs))
    med_spread = float(np.nanmedian(np.abs(meds - np.nanmedian(meds))))
 
    return {
        "tile_bg_sigma_ratio": ratio,
        "tile_bg_med_spread": med_spread,
        "global_bg_sigma_e_per_s": float(g_sig),
        "global_bg_med_e_per_s": float(g_med),
        "tiles_used": int(len(sigs)),
    }
 
 
def global_gradient(data: np.ndarray, wht: Optional[np.ndarray], step: int = 8) -> float:
    """
    Coarse global background gradient proxy:
    - downsample
    - compute median(|dx| + |dy|) across valid pixels
    """
    if wht is not None:
        valid = np.isfinite(data) & np.isfinite(wht) & (wht > 0)
    else:
        valid = np.isfinite(data)
 
    x = downsample2d(data, step)
    vm = downsample2d(valid.astype(np.uint8), step) > 0
    x = np.where(vm & np.isfinite(x), x, np.nan)
 
    dx = np.nanmedian(np.abs(np.diff(x, axis=1)))
    dy = np.nanmedian(np.abs(np.diff(x, axis=0)))
    if not np.isfinite(dx) and not np.isfinite(dy):
        return float("nan")
    return float(np.nan_to_num(dx, nan=0.0) + np.nan_to_num(dy, nan=0.0))
 
 
def sharpness_proxy(data: np.ndarray, wht: Optional[np.ndarray], step: int = 2) -> float:
    """
    High-frequency energy proxy for mosaic sharpness.
    """
    x = downsample2d(data, step)
    if wht is not None:
        v = downsample2d((np.isfinite(data) & np.isfinite(wht) & (wht > 0)).astype(np.uint8), step) > 0
    else:
        v = downsample2d(np.isfinite(data).astype(np.uint8), step) > 0
 
    x = np.where(v & np.isfinite(x), x, 0.0)
 
    lap = (
        -4.0 * x
        + np.roll(x, 1, 0) + np.roll(x, -1, 0)
        + np.roll(x, 1, 1) + np.roll(x, -1, 1)
    )
 
    hf = np.median(np.abs(lap[v])) if np.any(v) else float("nan")
    sig = robust_sigma(x[v]) if np.any(v) else float("nan")
    if not np.isfinite(hf) or not np.isfinite(sig) or sig <= 0:
        return float("nan")
    return float(hf / sig)
 
 
# ----------------------------
# WCS plausibility (lightweight)
# ----------------------------
 
def wcs_checks(model) -> Dict[str, float]:
    """
    Lightweight checks: ensure essential WCS meta exists; compute pixel scale from CD matrix if present.
    """
    out = {
        "wcs_present": 0.0,
        "pixel_scale_arcsec": float("nan"),
        "pixel_scale_frac_dev": float("nan"),
    }
 
    try:
        wcs = getattr(model.meta, "wcs", None)
        if wcs is None:
            return out
        out["wcs_present"] = 1.0
    except Exception:
        return out
 
    # Try to infer pixel scale from FITS WCS keywords if present
    try:
        # data models often store FITS-like WCS in meta.wcsinfo
        w = getattr(model.meta, "wcsinfo", None)
        if w is None:
            return out
 
        cd11 = safe_float(getattr(w, "cd1_1", None))
        cd12 = safe_float(getattr(w, "cd1_2", None))
        cd21 = safe_float(getattr(w, "cd2_1", None))
        cd22 = safe_float(getattr(w, "cd2_2", None))
 
        if None in (cd11, cd12, cd21, cd22):
            # try CDELT if CD missing
            cdelt1 = safe_float(getattr(w, "cdelt1", None))
            cdelt2 = safe_float(getattr(w, "cdelt2", None))
            if cdelt1 is None or cdelt2 is None:
                return out
            # degrees per pixel -> arcsec per pixel
            scale = 3600.0 * float((abs(cdelt1) + abs(cdelt2)) / 2.0)
        else:
            # scale in deg/pixel; approximate using sqrt(|det(CD)|)
            det = cd11 * cd22 - cd12 * cd21
            scale = 3600.0 * math.sqrt(abs(det))
 
        out["pixel_scale_arcsec"] = float(scale)
    except Exception:
        return out
 
    return out
 
 
# ----------------------------
# Scoring & classification
# ----------------------------
 
def compute_scores(metrics: Dict[str, float], thr: NircamMosaicThresholds) -> Dict[str, float]:
    """
    RTT-style scores:
    - noise: driven by global bg sigma (normalized loosely through seam contrast and tile ratios)
    - risk: seam + wcs issues + sharpness outliers + gradient
    - stability: inverse of risk/noise
    - internal: wcs presence + finite core metrics
    - memory: metadata completeness proxy (via internal)
    """
    seam_contrast = metrics.get("seam_contrast", float("nan"))
    seam_abs = metrics.get("seam_abs_e_per_s", float("nan"))
    tile_ratio = metrics.get("tile_bg_sigma_ratio", float("nan"))
    grad = metrics.get("global_gradient_e_per_s_per_pix", float("nan"))
    sharp = metrics.get("sharpness", float("nan"))
    wcs_present = metrics.get("wcs_present", 0.0)
    pix_scale_dev = metrics.get("pixel_scale_frac_dev", float("nan"))
 
    # Noise proxy: tile sigma non-uniformity + gradient
    noise = 0.0
    if np.isfinite(tile_ratio):
        noise += clamp((tile_ratio - 1.0) / max(1e-6, thr.max_tile_bg_sigma_ratio - 1.0))
    else:
        noise += 0.4
    if np.isfinite(grad):
        noise += 0.6 * clamp(grad / max(1e-6, thr.max_global_gradient_e_per_s_per_pix))
    else:
        noise += 0.3
    noise = clamp(noise)
 
    # Risk proxy: seam boundary crossing + sharpness out-of-family + missing WCS
    seam_term = 0.0
    if np.isfinite(seam_contrast):
        seam_term = clamp(seam_contrast / max(1e-6, thr.max_seam_contrast))
    elif np.isfinite(seam_abs):
        seam_term = clamp(seam_abs / max(1e-6, thr.max_seam_abs_e_per_s))
    else:
        seam_term = 0.5
 
    sharp_term = 0.0
    if np.isfinite(sharp):
        if sharp < thr.min_sharpness:
            sharp_term = clamp((thr.min_sharpness - sharp) / thr.min_sharpness)
        elif sharp > thr.max_sharpness:
            sharp_term = clamp((sharp - thr.max_sharpness) / thr.max_sharpness)
 
    wcs_term = 1.0 - float(wcs_present)
 
    pix_term = 0.0
    if np.isfinite(pix_scale_dev):
        pix_term = clamp(pix_scale_dev / max(1e-6, thr.max_pixel_scale_frac))
 
    risk = clamp(0.45 * seam_term + 0.25 * wcs_term + 0.20 * sharp_term + 0.10 * pix_term)
 
    stability = clamp(1.0 - (0.55 * risk + 0.45 * noise))
    internal = clamp(0.6 * float(wcs_present) + 0.4 * float(np.isfinite(seam_contrast) or np.isfinite(seam_abs)))
    memory = internal
 
    coherence = clamp(0.40 * stability + 0.25 * internal + 0.15 * memory - 0.10 * noise - 0.10 * risk)
 
    return {
        "noise": float(noise),
        "risk": float(risk),
        "stability": float(stability),
        "internal": float(internal),
        "memory": float(memory),
        "coherence": float(coherence),
    }
 
 
def zone_and_why(metrics: Dict[str, float], scores: Dict[str, float], thr: NircamMosaicThresholds) -> Tuple[str, str]:
    reasons = []
 
    seam_contrast = metrics.get("seam_contrast", float("nan"))
    seam_abs = metrics.get("seam_abs_e_per_s", float("nan"))
    tile_ratio = metrics.get("tile_bg_sigma_ratio", float("nan"))
    grad = metrics.get("global_gradient_e_per_s_per_pix", float("nan"))
    sharp = metrics.get("sharpness", float("nan"))
    wcs_present = metrics.get("wcs_present", 0.0)
 
    if np.isfinite(seam_contrast) and seam_contrast > thr.max_seam_contrast:
        reasons.append(f"Seam contrast {seam_contrast:.2f} > {thr.max_seam_contrast}")
    elif np.isfinite(seam_abs) and seam_abs > thr.max_seam_abs_e_per_s:
        reasons.append(f"Seam abs {seam_abs:.2f} > {thr.max_seam_abs_e_per_s}")
 
    if np.isfinite(tile_ratio) and tile_ratio > thr.max_tile_bg_sigma_ratio:
        reasons.append(f"Tile sigma ratio {tile_ratio:.2f} > {thr.max_tile_bg_sigma_ratio}")
 
    if np.isfinite(grad) and grad > thr.max_global_gradient_e_per_s_per_pix:
        reasons.append(f"Gradient {grad:.3f} > {thr.max_global_gradient_e_per_s_per_pix}")
 
    if np.isfinite(sharp) and (sharp < thr.min_sharpness or sharp > thr.max_sharpness):
        reasons.append(f"Sharpness {sharp:.3f} out of [{thr.min_sharpness},{thr.max_sharpness}]")
 
    if float(wcs_present) < 0.5:
        reasons.append("WCS missing")
 
    if scores["risk"] >= thr.verge_risk_min:
        why = "; ".join(reasons) if reasons else "Risk high: seam/uniformity boundary crossed"
        return "Transit Verge", why
 
    if scores["stability"] >= thr.calm_stability_min and scores["risk"] < 0.30:
        return "Lagrange Calm", "Stable field; seams and background within bounds"
 
    if scores["noise"] < 0.30 and scores["internal"] > 0.70:
        return "Deep Quiet", "Exceptionally coherent mosaic (low structure, strong metadata coherence)"
 
    why = "; ".join(reasons) if reasons else "Nominal mosaic; trend watch"
    return "Echo Belt", why
 
 
# ----------------------------
# Optional PNG overlay
# ----------------------------
 
def write_seam_png(png_path: Path, data: np.ndarray, seam_mask: np.ndarray, wht: Optional[np.ndarray]) -> None:
    """
    Write a quicklook PNG with seam mask overlay (no photometric claims).
    Requires matplotlib. Off by default.
    """
    try:
        import matplotlib.pyplot as plt
    except Exception as e:
        raise RuntimeError("matplotlib required for --png overlays") from e
 
    # Build a display image using robust stretch on background-ish pixels
    if wht is not None:
        valid = np.isfinite(data) & np.isfinite(wht) & (wht > 0)
    else:
        valid = np.isfinite(data)
 
    v = data[valid]
    if v.size < 1000:
        img = np.zeros_like(data, dtype=float)
        vmin, vmax = 0, 1
    else:
        p1, p99 = np.percentile(v, 1), np.percentile(v, 99)
        vmin, vmax = float(p1), float(p99)
        img = np.clip(data, vmin, vmax)
 
    # Downsample for speed and file size
    step = 4 if max(data.shape) > 3000 else 2
    img_s = downsample2d(img, step)
    seam_s = downsample2d(seam_mask.astype(np.uint8), step) > 0
 
    plt.figure(figsize=(10, 8))
    plt.imshow(img_s, cmap="gray", vmin=vmin, vmax=vmax)
    # seam overlay in red
    overlay = np.zeros((*img_s.shape, 4), dtype=float)
    overlay[seam_s, :] = [1.0, 0.1, 0.1, 0.65]
    plt.imshow(overlay)
    plt.title("RTT seam overlay (quicklook)")
    plt.axis("off")
    plt.tight_layout()
    plt.savefig(png_path, dpi=150)
    plt.close()
 
 
# ----------------------------
# Outputs
# ----------------------------
 
def write_json(out_path: Path, report: Dict[str, Any]) -> None:
    out_path.write_text(json.dumps(report, indent=2), encoding="utf-8")
 
 
def write_csv_row(csv_path: Path, row: Dict[str, Any]) -> None:
    """
    Writes/append a one-row CSV with a stable column order.
    """
    cols = list(row.keys())
    exists = csv_path.exists()
 
    with csv_path.open("a", newline="", encoding="utf-8") as f:
        w = csv.DictWriter(f, fieldnames=cols)
        if not exists:
            w.writeheader()
        w.writerow(row)
 
 
# ----------------------------
# Main
# ----------------------------
 
def main() -> None:
    ap = argparse.ArgumentParser(description="RTT sidecar QA for JWST NIRCam i2d mosaics (seam/coadd coherence).")
    ap.add_argument("fits", nargs="+", help="Input i2d.fits files (one or many).")
    ap.add_argument("--band", default=None, help="SW or LW. If omitted, tries to infer from filter.")
    ap.add_argument("--grid", type=int, default=8, help="Grid tiles per axis for segmentation fallback (default 8).")
    ap.add_argument("--seam-band", type=int, default=8, help="Seam evaluation band width in pixels (default 8).")
    ap.add_argument("--out", default=None, help="Output JSON path (single input only). Default: <input>.rtt.json")
    ap.add_argument("--csv-out", default=None, help="Optional CSV summary path (aggregates multiple inputs).")
    ap.add_argument("--png", action="store_true", help="Write seam overlay PNG quicklook (default off).")
    ap.add_argument("--png-dir", default=None, help="Directory for PNG outputs (default alongside input).")
    args = ap.parse_args()
 
    inputs = [Path(p) for p in args.fits]
 
    # Single-file explicit JSON output path only makes sense for one input
    if args.out and len(inputs) != 1:
        raise ValueError("--out only valid when processing a single FITS file")
 
    for in_path in inputs:
        with datamodels.open(str(in_path)) as model:
            data, wht = get_arrays(model)
 
            band = args.band
            if band is None:
                band = infer_band_from_meta(model)
            if band is None:
                raise ValueError(f"Could not infer band for {in_path.name}; specify --band SW or --band LW")
 
            thr = pick_thresholds(band)
 
            if wht is not None:
                valid = np.isfinite(data) & np.isfinite(wht) & (wht > 0)
            else:
                valid = np.isfinite(data)
 
            tile_id, seg_info = segmentation_from_weights(wht, grid=args.grid)
            seam_from_tiles = seam_mask_from_tile_ids(tile_id, valid)
 
            # If segmentation_from_weights included a weight-boundary mask idea, we still
            # use tile seams as primary seam loci for v0.1 (stable + explainable).
            seam_mask = seam_from_tiles
 
            seam = seam_metrics(data, wht, tile_id, seam_mask, band_width=args.seam_band)
            tiles = tile_stats(data, wht, tile_id)
            grad = global_gradient(data, wht, step=8)
            sharp = sharpness_proxy(data, wht, step=2)
            wcs = wcs_checks(model)
 
            # Pixel scale dev (if pixel scale exists, compare to nominal NIRCam SW/LW scale)
            # Nominal scales (approx): SW ~0.031 arcsec/pix, LW ~0.063 arcsec/pix
            nominal = 0.031 if band.upper() == "SW" else 0.063
            pix_scale = wcs.get("pixel_scale_arcsec", float("nan"))
            pix_dev = float(abs(pix_scale - nominal) / nominal) if np.isfinite(pix_scale) else float("nan")
            wcs["pixel_scale_frac_dev"] = pix_dev
 
            metrics: Dict[str, float] = {
                "band": band,
                "valid_frac": float(np.count_nonzero(valid) / valid.size),
                "seam_abs_e_per_s": seam.get("seam_abs_e_per_s", float("nan")),
                "seam_contrast": seam.get("seam_contrast", float("nan")),
                "seam_density": seam.get("seam_density", float("nan")),
                "bg_sigma_e_per_s": seam.get("bg_sigma_e_per_s", float("nan")),
                "tile_bg_sigma_ratio": tiles.get("tile_bg_sigma_ratio", float("nan")),
                "tile_bg_med_spread": tiles.get("tile_bg_med_spread", float("nan")),
                "global_bg_sigma_e_per_s": tiles.get("global_bg_sigma_e_per_s", float("nan")),
                "global_bg_med_e_per_s": tiles.get("global_bg_med_e_per_s", float("nan")),
                "global_gradient_e_per_s_per_pix": float(grad),
                "sharpness": float(sharp),
                "wcs_present": float(wcs.get("wcs_present", 0.0)),
                "pixel_scale_arcsec": float(wcs.get("pixel_scale_arcsec", float("nan"))),
                "pixel_scale_frac_dev": float(wcs.get("pixel_scale_frac_dev", float("nan"))),
                "tiles_used": float(tiles.get("tiles_used", float("nan"))),
            }
 
            scores = compute_scores(metrics, thr)
            zone, why = zone_and_why(metrics, scores, thr)
 
            # Lightweight metadata capture (non-sensitive)
            meta: Dict[str, Any] = {}
            try:
                meta = {
                    "instrument": getattr(model.meta.instrument, "name", None),
                    "detector": getattr(model.meta.instrument, "detector", None),
                    "filter": getattr(model.meta.instrument, "filter", None),
                    "pupil": getattr(model.meta.instrument, "pupil", None),
                    "program": getattr(model.meta.observation, "program_number", None),
                    "obs": getattr(model.meta.observation, "observation_number", None),
                    "visit": getattr(model.meta.observation, "visit_number", None),
                    "exp_type": getattr(model.meta.exposure, "type", None),
                }
            except Exception:
                pass
 
            report: Dict[str, Any] = {
                "rtt_version": "0.1-nircam-i2d",
                "input": str(in_path),
                "zone": zone,
                "why": why,
                "metrics": metrics,
                "scores": scores,
                "thresholds": asdict(thr),
                "segmentation": seg_info,
                "meta": meta,
                "notes": "Sidecar QA only; does not modify calibrated arrays. Mosaic-level seam/coadd coherence lens.",
            }
 
        # Write JSON
        out_json = Path(args.out) if args.out else in_path.with_suffix(in_path.suffix + ".rtt.json")
        write_json(out_json, report)
 
        # Write per-file CSV summary unless batch CSV provided
        row = {
            "input": str(in_path),
            "band": metrics["band"],
            "zone": zone,
            "why": why,
            "coherence": scores["coherence"],
            "risk": scores["risk"],
            "noise": scores["noise"],
            "stability": scores["stability"],
            "seam_contrast": metrics["seam_contrast"],
            "seam_abs_e_per_s": metrics["seam_abs_e_per_s"],
            "tile_bg_sigma_ratio": metrics["tile_bg_sigma_ratio"],
            "global_gradient_e_per_s_per_pix": metrics["global_gradient_e_per_s_per_pix"],
            "sharpness": metrics["sharpness"],
            "wcs_present": metrics["wcs_present"],
            "pixel_scale_arcsec": metrics["pixel_scale_arcsec"],
            "pixel_scale_frac_dev": metrics["pixel_scale_frac_dev"],
        }
 
        if args.csv_out:
            write_csv_row(Path(args.csv_out), row)
        else:
            write_csv_row(in_path.with_suffix(in_path.suffix + ".rtt.csv"), row)
 
        # Optional PNG
        if args.png:
            png_dir = Path(args.png_dir) if args.png_dir else in_path.parent
            png_dir.mkdir(parents=True, exist_ok=True)
            png_path = png_dir / (in_path.name + ".rtt_seams.png")
 
            # Re-open for arrays (avoid keeping model open for plotting)
            with datamodels.open(str(in_path)) as model:
                data, wht = get_arrays(model)
                if wht is not None:
                    valid = np.isfinite(data) & np.isfinite(wht) & (wht > 0)
                else:
                    valid = np.isfinite(data)
                tile_id, _ = segmentation_from_weights(wht, grid=args.grid)
                seam_mask = seam_mask_from_tile_ids(tile_id, valid)
                write_seam_png(png_path, data, seam_mask, wht)
 
        print(f"[RTT] {in_path.name}: zone={zone} coherence={scores['coherence']:.2f} -> {out_json}")
 
 
if __name__ == "__main__":
    main()

scores = compute_scores(metrics, thr)
zone, why = zone_and_why(metrics, scores, thr)