"""Binary 2-D metrics with explicit empty-mask and pixel-center conventions.

Distances use all inner boundary pixel centers, four-neighbor erosion with
outside background. HD95 is max(directional 95th percentiles), method='linear'.
Pairwise distances are a small-mask teaching reference, not a volume backend.
"""
from __future__ import annotations

import cv2
import numpy as np
from numpy.typing import NDArray

Mask = NDArray[np.bool_]
FloatArray = NDArray[np.float64]


def validate_pair(pred: Mask, target: Mask) -> None:
    """Require matching, nonempty 2-D boolean grids; foreground may be empty."""
    if pred.dtype != np.bool_ or target.dtype != np.bool_:
        raise ValueError("Masks must have bool dtype")
    if pred.ndim != 2 or pred.shape != target.shape or pred.size == 0:
        raise ValueError("Masks must be matching nonempty 2-D grids")

# region overlap
def overlap(pred: Mask, target: Mask) -> tuple[float, float]:
    """Return IoU, Dice; both empty=1, exactly one empty=0."""
    validate_pair(pred, target)
    intersection: int = int(np.count_nonzero(pred & target))
    union: int = int(np.count_nonzero(pred | target))
    total: int = int(np.count_nonzero(pred)) + int(np.count_nonzero(target))
    if union == 0:
        return 1.0, 1.0
    iou: float = intersection / union
    dice: float = 2.0 * intersection / total
    return iou, dice
# endregion overlap

# region boundary
def boundary(mask: Mask) -> Mask:
    """Return inner boundary pixels using a center+four-neighbor cross."""
    validate_pair(mask, mask)
    cross: NDArray[np.uint8] = np.array(
        [[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=np.uint8,
    )
    interior: Mask = cv2.erode(
        mask.astype(np.uint8), cross,
        borderType=cv2.BORDER_CONSTANT, borderValue=0,
    ) != 0
    edge: Mask = mask & ~interior
    return edge
# endregion boundary

# region distances
def directed_distances(
    source: Mask, destination: Mask, spacing: tuple[float, float] = (1.0, 1.0),
) -> FloatArray:
    """Each source boundary center -> nearest destination center distance.

    spacing=(row, column), in pixels or a declared physical unit. This dense
    reference rejects products above 2 million pairs before allocating them.
    Caller handles empty boundaries; output float64 has shape (source_count,).
    """
    validate_pair(source, destination)
    scale: FloatArray = np.asarray(spacing, dtype=np.float64)
    if scale.shape != (2,) or not np.isfinite(scale).all() or np.any(scale <= 0):
        raise ValueError("Need positive finite row/column spacing")
    a: FloatArray = np.argwhere(boundary(source)).astype(np.float64) * scale
    b: FloatArray = np.argwhere(boundary(destination)).astype(np.float64) * scale
    if len(a) == 0 or len(b) == 0:
        raise ValueError("Directed distances need two nonempty boundaries")
    if len(a) * len(b) > 2_000_000:
        raise ValueError("Use a distance-transform backend for large masks")
    delta: FloatArray = a[:, None, :] - b[None, :, :]
    squared: FloatArray = np.sum(delta * delta, axis=-1)
    nearest: FloatArray = np.sqrt(squared.min(axis=1))
    return nearest
# endregion distances

# region hausdorff
def hausdorff(
    pred: Mask, target: Mask, spacing: tuple[float, float] = (1.0, 1.0),
) -> tuple[float, float]:
    """Return symmetric HD and max-directional HD95; empty=0 or +inf."""
    validate_pair(pred, target)
    if len(spacing) != 2 or not all(np.isfinite(s) and s > 0 for s in spacing):
        raise ValueError("Need positive finite row/column spacing")
    if not pred.any() and not target.any():
        return 0.0, 0.0
    if not pred.any() or not target.any():
        return float("inf"), float("inf")
    forward: FloatArray = directed_distances(pred, target, spacing)
    backward: FloatArray = directed_distances(target, pred, spacing)
    hd: float = float(max(forward.max(), backward.max()))
    hd95: float = float(max(
        np.percentile(forward, 95, method="linear"),
        np.percentile(backward, 95, method="linear"),
    ))
    return hd, hd95
# endregion hausdorff


def metric_cases() -> dict[str, tuple[Mask, Mask]]:
    """All pairs are (prediction, target); fixed 12x12 grids, no randomness."""
    target: Mask = np.zeros((12, 12), dtype=np.bool_)
    target[3:9, 3:9] = True
    shifted: Mask = np.zeros_like(target)
    shifted[3:9, 4:10] = True
    outlier: Mask = target.copy()
    outlier[0, 0] = True
    disjoint: Mask = np.zeros_like(target)
    disjoint[:2, :2] = True
    empty: Mask = np.zeros_like(target)
    return {
        "exact": (target.copy(), target), "shifted": (shifted, target),
        "outlier": (outlier, target), "disjoint": (disjoint, target),
        "both-empty": (empty.copy(), empty), "one-empty": (empty, target),
    }