"""Polygon-based morphometric traits for a single segmented tuber.

The width and height used by downstream stages are the extents of the
polygon bounding box. The aspect_ratio field is bbox max/min and is not
used as the published axis-ratio; a PCA-based ratio is derived from
eccentricity in the modeling stage.
"""

from __future__ import annotations

import numpy as np


def close_polygon(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    if x.size == 0 or y.size == 0:
        return x, y
    if x[0] != x[-1] or y[0] != y[-1]:
        x = np.concatenate([x, x[:1]])
        y = np.concatenate([y, y[:1]])
    return x, y


def polygon_area_px(x: np.ndarray, y: np.ndarray) -> float:
    """Shoelace formula. Eq. (2) in the paper."""
    x, y = close_polygon(x, y)
    if x.size < 3:
        return float("nan")
    area = 0.0
    for i in range(len(x) - 1):
        area += x[i] * y[i + 1] - x[i + 1] * y[i]
    return abs(area) / 2.0


def polygon_perimeter_px(x: np.ndarray, y: np.ndarray) -> float:
    """Sum of Euclidean distances between consecutive vertices. Eq. (1) in the paper."""
    x, y = close_polygon(x, y)
    if x.size < 2:
        return float("nan")
    dx = np.diff(x)
    dy = np.diff(y)
    return float(np.sqrt(dx**2 + dy**2).sum())


def polygon_circularity(x: np.ndarray, y: np.ndarray) -> float:
    """4·π·A / P²."""
    A = polygon_area_px(x, y)
    P = polygon_perimeter_px(x, y)
    if not np.isfinite(A) or not np.isfinite(P) or P <= 0:
        return float("nan")
    return float(4.0 * np.pi * A / (P**2))


def polygon_eccentricity(x: np.ndarray, y: np.ndarray) -> float:
    """sqrt(1 - λ_min / λ_max) of the polygon-vertex covariance matrix."""
    if len(x) < 3:
        return float("nan")
    coords = np.column_stack([x, y])
    coords_centered = coords - coords.mean(axis=0, keepdims=True)
    cov_mat = np.cov(coords_centered, rowvar=False)
    eigvals = np.linalg.eigvalsh(cov_mat)
    eigvals = np.sort(eigvals)[::-1]
    if eigvals.size < 2 or eigvals[0] <= 0:
        return float("nan")
    return float(np.sqrt(1.0 - eigvals[1] / eigvals[0]))


def aspect_ratio_bbox(width: float, height: float) -> float:
    if width <= 0 or height <= 0:
        return float("nan")
    return float(max(width, height) / min(width, height))


def roundness_bbox(width: float, height: float) -> float:
    if width <= 0 or height <= 0:
        return float("nan")
    area_bbox = width * height
    max_dim = max(width, height)
    return float(4.0 * area_bbox / (np.pi * max_dim**2))


def polygon_shape_metrics(x: np.ndarray, y: np.ndarray) -> dict[str, float]:
    """All per-tuber pixel-space metrics in a single call.

    Returns: area_poly_px, perimeter_px, circularity, eccentricity,
    width_px, height_px, aspect_ratio, roundness_bbox.
    """
    area_px = polygon_area_px(x, y)
    perim_px = polygon_perimeter_px(x, y)
    circ = polygon_circularity(x, y)
    ecc = polygon_eccentricity(x, y)

    width_px = float(x.max() - x.min()) if x.size > 0 else float("nan")
    height_px = float(y.max() - y.min()) if y.size > 0 else float("nan")
    ar = aspect_ratio_bbox(width_px, height_px)
    rb = roundness_bbox(width_px, height_px)

    return {
        "area_poly_px": area_px,
        "perimeter_px": perim_px,
        "circularity": circ,
        "eccentricity": ecc,
        "width_px": width_px,
        "height_px": height_px,
        "aspect_ratio": ar,
        "roundness_bbox": rb,
    }