"""Validation metrics — segmentation and trajectory accuracy.

Evidence a surgeon / IRB will ask for before clinical validation: how close are
the segmentations to ground truth (Dice, Jaccard, Hausdorff-95, average surface
distance, volume error), and how close are planned trajectories to a reference
(entry/target error in mm, angular error). Pure NumPy + scipy; unit-tested.

Use with `scripts/validate_segmentation.py` to score a registered model against
annotated masks and emit a metrics report.

Research software only. Not a medical device.
"""
from __future__ import annotations
import numpy as np


def _bool(m):
    return np.asarray(getattr(m, "data", m)) > 0


def _spacing(m, spacing=None):
    if spacing is not None:
        return np.asarray(spacing, float)
    sp = getattr(m, "spacing", None)
    return np.asarray(sp, float) if sp is not None else np.ones(3)


def dice(pred, truth) -> float:
    """Dice similarity coefficient (2|A∩B| / (|A|+|B|)). 1.0 if both empty."""
    a, b = _bool(pred), _bool(truth)
    s = int(a.sum()) + int(b.sum())
    return 1.0 if s == 0 else 2.0 * int((a & b).sum()) / s


def jaccard(pred, truth) -> float:
    a, b = _bool(pred), _bool(truth)
    u = int((a | b).sum())
    return 1.0 if u == 0 else int((a & b).sum()) / u


def volume_cc(mask, spacing=None) -> float:
    return round(float(_bool(mask).sum()) * float(np.prod(_spacing(mask, spacing))) / 1000.0, 2)


def _surface(m):
    from scipy import ndimage
    return m & ~ndimage.binary_erosion(m)


def surface_distances(pred, truth, spacing=None):
    """Symmetric surface distances (mm): (pred→truth, truth→pred) arrays."""
    from scipy import ndimage
    a, b = _bool(pred), _bool(truth)
    sp = _spacing(pred, spacing)
    if a.sum() == 0 or b.sum() == 0:
        return None
    sa, sb = _surface(a), _surface(b)
    dtb = ndimage.distance_transform_edt(~sb, sampling=sp)
    dta = ndimage.distance_transform_edt(~sa, sampling=sp)
    return dtb[sa], dta[sb]


def hausdorff95(pred, truth, spacing=None) -> float:
    r = surface_distances(pred, truth, spacing)
    if r is None:
        return float("inf")
    return float(np.percentile(np.concatenate(r), 95))


def average_surface_distance(pred, truth, spacing=None) -> float:
    r = surface_distances(pred, truth, spacing)
    if r is None:
        return float("inf")
    da, db = r
    return float((da.mean() + db.mean()) / 2.0)


def segmentation_metrics(pred, truth, spacing=None) -> dict:
    """Full segmentation report vs ground truth."""
    sp = _spacing(pred, spacing)
    vp, vt = volume_cc(pred, sp), volume_cc(truth, sp)
    return {
        "dice": round(dice(pred, truth), 4),
        "jaccard": round(jaccard(pred, truth), 4),
        "volume_pred_cc": vp,
        "volume_truth_cc": vt,
        "volume_error_cc": round(vp - vt, 2),
        "abs_volume_error_pct": round(100 * abs(vp - vt) / vt, 1) if vt else float("inf"),
        "hausdorff95_mm": round(hausdorff95(pred, truth, sp), 2),
        "assd_mm": round(average_surface_distance(pred, truth, sp), 2),
    }


def trajectory_error(pred, ref) -> dict:
    """Entry/target error (mm) and angular error (deg) vs a reference trajectory.
    Accepts objects with .entry/.target or (entry, target) pairs."""
    def ee(t):
        return (np.asarray(t.entry, float), np.asarray(t.target, float)) if hasattr(t, "entry") \
            else (np.asarray(t[0], float), np.asarray(t[1], float))
    pe, pt = ee(pred)
    re, rt = ee(ref)
    du = pt - pe; dv = rt - re
    nu, nv = np.linalg.norm(du), np.linalg.norm(dv)
    ang = 0.0
    if nu > 0 and nv > 0:
        ang = float(np.degrees(np.arccos(np.clip(np.dot(du, dv) / (nu * nv), -1, 1))))
    return {"entry_error_mm": round(float(np.linalg.norm(pe - re)), 2),
            "target_error_mm": round(float(np.linalg.norm(pt - rt)), 2),
            "angle_error_deg": round(ang, 2)}


def summarize(metrics_list, keys=("dice", "hausdorff95_mm", "assd_mm", "abs_volume_error_pct")) -> dict:
    """Mean ± std over a list of per-case metric dicts (for a dataset report)."""
    out = {"n": len(metrics_list)}
    for k in keys:
        vals = [m[k] for m in metrics_list if k in m and np.isfinite(m[k])]
        if vals:
            out[k] = {"mean": round(float(np.mean(vals)), 3), "std": round(float(np.std(vals)), 3),
                      "min": round(float(np.min(vals)), 3), "max": round(float(np.max(vals)), 3)}
    return out
