"""輝度分布「形」評価指標モジュール．

現状の均一性指標（CoV / 勾配 / C/P / 動径 min/max）は不均一性の絶対量しか測れず，
Y 軸単調勾配と中央ピーク楕円のように「形が構造的に異なる」パターンを区別できない．
本モジュールは「形」を直接評価する 5 つの純関数を提供する．

References:
    [1] TECH_11_形状指標の見直しと新指標導入
        docs/05_TECH/TECH_11_形状指標の見直しと新指標導入.md
"""

from __future__ import annotations

import numpy


def calc_intensity_centroid(
    image: numpy.ndarray,
    mask: numpy.ndarray | None = None,
) -> dict[str, float]:
    """輝度重心を計算する．

    各ピクセルの輝度値を重みとして ROI 内の重心座標を算出する．
    重心が ROI 中心から大きくずれていれば，均一でない輝度分布や
    非対称な光源配置を示す．

    Args:
        image: 2D (H, W) 輝度マップ．負値は ReLU で 0 に．
        mask: 2D bool マスク（image と同 shape）．True のピクセルのみ計算対象．
            None なら全画素を対象とする．

    Returns:
        以下のキーを持つ辞書:
        - ``"centroid_x_px"`` (float): ROI 座標系での重心 X（列方向，左端 = 0）．
        - ``"centroid_y_px"`` (float): ROI 座標系での重心 Y（行方向，上端 = 0）．
        - ``"centroid_x_norm"`` (float): ROI 中心を原点とした正規化 X（-1〜+1）．
        - ``"centroid_y_norm"`` (float): ROI 中心を原点とした正規化 Y（-1〜+1）．
        - ``"centroid_offset_norm"`` (float): 正規化座標での中心からの距離 √(x²+y²)．

    Raises:
        ValueError: mask の shape が image と異なる場合．
        ValueError: 有効ピクセルの輝度総和が 0 の場合（重心計算不能）．
    """
    if image.ndim != 2:
        msg = f"image は 2D 配列を期待します（ndim={image.ndim}）"
        raise ValueError(msg)

    if mask is not None:
        if mask.shape != image.shape:
            msg = (
                f"mask の shape {mask.shape} が image の shape {image.shape} と異なります"
            )
            raise ValueError(msg)

    h, w = image.shape

    # 負値は 0 に（ReLU）
    img = numpy.maximum(image.astype(numpy.float64), 0.0)

    if mask is not None:
        img = numpy.where(mask, img, 0.0)

    total_weight = float(img.sum())
    if total_weight == 0.0:
        msg = "有効ピクセルの輝度総和が 0 のため重心を計算できません"
        raise ValueError(msg)

    # グリッド座標（ピクセル中心）
    ys = numpy.arange(h, dtype=numpy.float64).reshape(-1, 1)
    xs = numpy.arange(w, dtype=numpy.float64).reshape(1, -1)

    centroid_x_px = float((img * xs).sum() / total_weight)
    centroid_y_px = float((img * ys).sum() / total_weight)

    # ROI 中心を原点とした正規化座標
    cx = (w - 1) / 2.0
    cy = (h - 1) / 2.0
    # 正規化: ±1 が ROI の端になるよう cx/cy で割る
    centroid_x_norm = (centroid_x_px - cx) / (cx if cx > 0.0 else 1.0)
    centroid_y_norm = (centroid_y_px - cy) / (cy if cy > 0.0 else 1.0)
    centroid_offset_norm = float(numpy.sqrt(centroid_x_norm**2 + centroid_y_norm**2))

    return {
        "centroid_x_px": centroid_x_px,
        "centroid_y_px": centroid_y_px,
        "centroid_x_norm": centroid_x_norm,
        "centroid_y_norm": centroid_y_norm,
        "centroid_offset_norm": centroid_offset_norm,
    }


def calc_angular_radial_decomposition(
    image: numpy.ndarray,
    n_radial_bins: int = 20,
    n_angular_bins: int = 36,
    mask: numpy.ndarray | None = None,
) -> dict[str, float]:
    """輝度マップを極座標で角度/半径方向に分解する．

    画像を画像中心（ROI 中心）を原点とする極座標 (r, θ) に変換し，
    放射方向の分散と角度方向の分散を分離することで「形」の方向性を数値化する．

    アルゴリズム:
        1. 画像中心を原点とする極座標に変換
        2. (r_bin, θ_bin) のグリッドで輝度を分類 → I[r_bin, θ_bin]
        3. angular_variance = mean_r[ Var_θ[I | r] ]  各 r での角度方向 var の半径平均
        4. radial_variance  = Var_r[ mean_θ[I | r] ]  角度平均した放射プロファイルの var
        5. angular_ratio = angular_variance / (angular_variance + radial_variance)

    物理的意味:
        - ``angular_ratio ≈ 0``: 完全同心円（各 r での角度変化が小 → 中央ピーク楕円に近い）
        - ``angular_ratio ≈ 0.5``: 方向性ありの典型的非均一パターン
        - ``angular_ratio ≈ 1``: 完全角度依存（強い方向性ムラ）
        - 中央ピーク楕円: angular_ratio が小さい（0.1〜0.3）
        - Y 単調勾配: angular_ratio が中程度（0.3〜0.5）

    Notes:
        TECH_11 で示した「Y 単調勾配 → angular_ratio 0.3〜0.5」は
        実シミュ画像（ノイズ込み・有限解像度）での典型値である．
        純粋な Y 軸方向均一勾配（横縞画像）を合成すると，極座標のビン集計上
        角度方向変動が支配的になるため angular_ratio は 0.5 を超える場合がある．
        「Y 軸勾配は同心円より angular_ratio が高い」という方向性は
        理論・合成画像ともに成立する．

    Args:
        image: 2D (H, W) 輝度マップ．負値は 0 として扱う．
        n_radial_bins: 半径方向のビン数．デフォルト 20．
        n_angular_bins: 角度方向のビン数（360° を等分）．デフォルト 36（各 10°）．
        mask: 2D bool マスク（image と同 shape）．True のピクセルのみ計算対象．
            None なら全画素を対象とする．

    Returns:
        以下のキーを持つ辞書:
        - ``"angular_variance"`` (float): 角度方向分散の平均
        - ``"radial_variance"`` (float): 放射方向分散
        - ``"angular_ratio"`` (float): 0〜1 の形指標（0=同心円，1=角度依存）．
          完全均一画像（全 variance = 0）では NaN を返す．

    Raises:
        ValueError: mask の shape が image と異なる場合．
    """
    if image.ndim != 2:
        msg = f"image は 2D 配列を期待します（ndim={image.ndim}）"
        raise ValueError(msg)

    if mask is not None:
        if mask.shape != image.shape:
            msg = (
                f"mask の shape {mask.shape} が image の shape {image.shape} と異なります"
            )
            raise ValueError(msg)

    h, w = image.shape
    img = numpy.maximum(image.astype(numpy.float64), 0.0)

    if mask is not None:
        img = numpy.where(mask, img, numpy.nan)

    # 画像中心を原点とした座標グリッド
    cy = (h - 1) / 2.0
    cx = (w - 1) / 2.0
    ys = numpy.arange(h, dtype=numpy.float64).reshape(-1, 1) - cy
    xs = numpy.arange(w, dtype=numpy.float64).reshape(1, -1) - cx

    # 極座標変換
    r = numpy.sqrt(xs**2 + ys**2)
    theta = numpy.arctan2(ys, xs)  # [-π, +π]

    # r を [0, r_max] の n_radial_bins ビンに割り当て
    r_max = float(r.max())
    if r_max == 0.0:
        return {
            "angular_variance": 0.0,
            "radial_variance": 0.0,
            "angular_ratio": float("nan"),
        }

    r_norm = r / r_max  # [0, 1]

    # θ を [-π, π] → [0, 2π] へシフトして [0, n_angular_bins) に bin 化
    theta_shifted = theta + numpy.pi  # [0, 2π]
    r_bin_indices = numpy.minimum(
        (r_norm * n_radial_bins).astype(int), n_radial_bins - 1
    )
    theta_bin_indices = numpy.minimum(
        (theta_shifted / (2.0 * numpy.pi) * n_angular_bins).astype(int),
        n_angular_bins - 1,
    )

    # 各 (r_bin, θ_bin) の輝度平均を集計
    # grid[r_bin, theta_bin] = 輝度平均（空セルは nan）
    grid = numpy.full((n_radial_bins, n_angular_bins), numpy.nan, dtype=numpy.float64)
    r_flat = r_bin_indices.ravel()
    t_flat = theta_bin_indices.ravel()
    v_flat = img.ravel()

    # 有効なピクセルだけを集計する
    valid_mask_flat = ~numpy.isnan(v_flat)
    valid_r = r_flat[valid_mask_flat]
    valid_t = t_flat[valid_mask_flat]
    valid_v = v_flat[valid_mask_flat]

    # numpy.add.at で各セルに値を累積し，カウントで割って平均を得る
    sum_grid = numpy.zeros((n_radial_bins, n_angular_bins), dtype=numpy.float64)
    count_grid = numpy.zeros((n_radial_bins, n_angular_bins), dtype=numpy.int64)
    numpy.add.at(sum_grid, (valid_r, valid_t), valid_v)
    numpy.add.at(count_grid, (valid_r, valid_t), 1)

    valid_cells = count_grid > 0
    grid[valid_cells] = sum_grid[valid_cells] / count_grid[valid_cells]

    # angular_variance: 各 r_bin について角度方向の分散を算出し，r 方向に平均
    angular_var_per_r = numpy.full(n_radial_bins, numpy.nan, dtype=numpy.float64)
    for r_idx in range(n_radial_bins):
        row = grid[r_idx, :]
        valid_vals = row[~numpy.isnan(row)]
        if valid_vals.size >= 2:
            angular_var_per_r[r_idx] = float(numpy.var(valid_vals, ddof=0))
        elif valid_vals.size == 1:
            angular_var_per_r[r_idx] = 0.0
        # else: 有効セルなし → nan のまま

    valid_ang_var = angular_var_per_r[~numpy.isnan(angular_var_per_r)]
    angular_variance = float(numpy.mean(valid_ang_var)) if valid_ang_var.size > 0 else 0.0

    # radial_variance: 各 r_bin について角度方向の平均を算出し，その r 方向の分散
    radial_mean_per_r = numpy.full(n_radial_bins, numpy.nan, dtype=numpy.float64)
    for r_idx in range(n_radial_bins):
        row = grid[r_idx, :]
        valid_vals = row[~numpy.isnan(row)]
        if valid_vals.size > 0:
            radial_mean_per_r[r_idx] = float(numpy.mean(valid_vals))

    valid_radial_means = radial_mean_per_r[~numpy.isnan(radial_mean_per_r)]
    radial_variance = (
        float(numpy.var(valid_radial_means, ddof=0)) if valid_radial_means.size >= 2 else 0.0
    )

    total_variance = angular_variance + radial_variance
    if total_variance == 0.0:
        angular_ratio = float("nan")
    else:
        angular_ratio = angular_variance / total_variance

    return {
        "angular_variance": angular_variance,
        "radial_variance": radial_variance,
        "angular_ratio": angular_ratio,
    }


def calc_ssim(
    image: numpy.ndarray,
    reference: numpy.ndarray,
    mask: numpy.ndarray | None = None,
    data_range: float | None = None,
) -> float:
    """構造的類似度 (Structural Similarity Index) を算出する．

    ``skimage.metrics.structural_similarity`` を使用する．
    形（空間構造）の類似性を 0〜1 で数値化する．1.0 で完全一致．

    Args:
        image: 評価対象 2D 輝度マップ（H, W）．
        reference: 参照画像 2D 輝度マップ（image と同 shape）．
        mask: 2D bool マスク．True ピクセルのみを計算対象とするため，
            mask 外を 0 で埋めた上で full=True モードの SSIM 全体値を用いる．
        data_range: 輝度値の最大レンジ（``image.max() - reference.max()`` の max）．
            None のとき ``max(image.max(), reference.max())`` を使う．

    Returns:
        SSIM 値（float, 0〜1）．

    Raises:
        ValueError: image と reference の shape が異なる場合．
        ValueError: mask の shape が image と異なる場合．
        ImportError: scikit-image がインストールされていない場合．
    """
    try:
        from skimage.metrics import structural_similarity  # noqa: PLC0415
    except ImportError as exc:
        msg = (
            "calc_ssim には scikit-image が必要です．"
            "pip install scikit-image>=0.21.0 でインストールしてください．"
        )
        raise ImportError(msg) from exc

    if image.ndim != 2 or reference.ndim != 2:
        msg = "image と reference は 2D 配列を期待します"
        raise ValueError(msg)

    if image.shape != reference.shape:
        msg = (
            f"image の shape {image.shape} と reference の shape {reference.shape} が異なります．"
            "事前に resize してください．"
        )
        raise ValueError(msg)

    if mask is not None:
        if mask.shape != image.shape:
            msg = (
                f"mask の shape {mask.shape} が image の shape {image.shape} と異なります"
            )
            raise ValueError(msg)

    img = image.astype(numpy.float64)
    ref = reference.astype(numpy.float64)

    if mask is not None:
        img = numpy.where(mask, img, 0.0)
        ref = numpy.where(mask, ref, 0.0)

    if data_range is None:
        data_range = float(max(img.max(), ref.max()))
        if data_range == 0.0:
            data_range = 1.0  # 全ゼロ画像でのゼロ除算を防ぐ

    ssim_val = structural_similarity(img, ref, data_range=data_range)
    return float(ssim_val)


def calc_ncc(
    image: numpy.ndarray,
    reference: numpy.ndarray,
    mask: numpy.ndarray | None = None,
) -> float:
    """正規化相関係数 (Normalized Cross-Correlation) を算出する．

    NCC = Σ((I - mean(I)) * (R - mean(R))) / sqrt(Σ(I - mean(I))² * Σ(R - mean(R))²)

    絶対輝度スケールに依存しない形の類似性指標．1.0 で完全相関（形が一致），
    0 で無相関，-1 で逆相関．SSIM より単純だが照明変動に対しロバスト．

    Args:
        image: 評価対象 2D 輝度マップ（H, W）．
        reference: 参照画像 2D 輝度マップ（image と同 shape）．
        mask: 2D bool マスク（image と同 shape）．True のピクセルのみ計算対象．
            None なら全画素を対象とする．

    Returns:
        NCC 値（float, -1〜+1）．
        分散がゼロの場合（均一画像）は NaN を返す．

    Raises:
        ValueError: image と reference の shape が異なる場合．
        ValueError: mask の shape が image と異なる場合．
    """
    if image.ndim != 2 or reference.ndim != 2:
        msg = "image と reference は 2D 配列を期待します"
        raise ValueError(msg)

    if image.shape != reference.shape:
        msg = (
            f"image の shape {image.shape} と reference の shape {reference.shape} が異なります"
        )
        raise ValueError(msg)

    if mask is not None:
        if mask.shape != image.shape:
            msg = (
                f"mask の shape {mask.shape} が image の shape {image.shape} と異なります"
            )
            raise ValueError(msg)

    img = image.astype(numpy.float64).ravel()
    ref = reference.astype(numpy.float64).ravel()

    if mask is not None:
        mask_flat = mask.ravel()
        img = img[mask_flat]
        ref = ref[mask_flat]

    if img.size == 0:
        return float("nan")

    img_centered = img - img.mean()
    ref_centered = ref - ref.mean()

    denom = float(
        numpy.sqrt(numpy.sum(img_centered**2) * numpy.sum(ref_centered**2))
    )
    if denom == 0.0:
        return float("nan")

    ncc_val = float(numpy.sum(img_centered * ref_centered) / denom)
    # 数値誤差で [-1, 1] の範囲をわずかに超える可能性があるためクリップ
    return float(numpy.clip(ncc_val, -1.0, 1.0))


def calc_intensity_pca(
    image: numpy.ndarray,
    mask: numpy.ndarray | None = None,
) -> dict[str, float]:
    """輝度分布の主成分分析で主軸方向と楕円率を計算する．

    各ピクセルを (x, y, weight=I(x,y)) とした weighted PCA で
    輝度分布の「伸び方向」と「偏心度」を数値化する．

    アルゴリズム:
        1. 輝度を重みとして重心座標 (cx, cy) を算出
        2. 重心からのオフセット (Δx, Δy) に輝度重みを掛けて共分散行列を構築
        3. 固有値分解 → 主軸方向（第 1 固有ベクトルの角度）と固有値

    物理的意味:
        - ``primary_axis_angle_deg ≈ 90°``: 主軸が Y 軸方向（Y 単調勾配の典型）
        - ``eccentricity ≈ 0``: 完全円形分布（等方的中央ピーク）
        - ``eccentricity ≈ 1``: 極端な楕円（強い方向性の輝度分布）

    Args:
        image: 2D (H, W) 輝度マップ．負値は 0 として扱う．
        mask: 2D bool マスク（image と同 shape）．True のピクセルのみ計算対象．
            None なら全画素を対象とする．

    Returns:
        以下のキーを持つ辞書:
        - ``"primary_axis_angle_deg"`` (float): 主軸の角度 [deg]（画像 x 軸基準，-90〜+90）．
        - ``"primary_eigenvalue"`` (float): 第 1 固有値（主軸方向の分散の大きさ）．
        - ``"secondary_eigenvalue"`` (float): 第 2 固有値（副軸方向の分散の大きさ）．
        - ``"eccentricity"`` (float): √(1 - secondary/primary)，0=完全円，1=極端な楕円．
          primary_eigenvalue が 0 の場合は 0 を返す．

    Raises:
        ValueError: mask の shape が image と異なる場合．
        ValueError: 有効ピクセルの輝度総和が 0 の場合（PCA 計算不能）．
    """
    if image.ndim != 2:
        msg = f"image は 2D 配列を期待します（ndim={image.ndim}）"
        raise ValueError(msg)

    if mask is not None:
        if mask.shape != image.shape:
            msg = (
                f"mask の shape {mask.shape} が image の shape {image.shape} と異なります"
            )
            raise ValueError(msg)

    h, w = image.shape
    img = numpy.maximum(image.astype(numpy.float64), 0.0)

    if mask is not None:
        img = numpy.where(mask, img, 0.0)

    # グリッド座標（ピクセル中心）
    ys = numpy.arange(h, dtype=numpy.float64).reshape(-1, 1)
    xs = numpy.arange(w, dtype=numpy.float64).reshape(1, -1)

    total_weight = float(img.sum())
    if total_weight == 0.0:
        msg = "有効ピクセルの輝度総和が 0 のため PCA を計算できません"
        raise ValueError(msg)

    # 重心
    cx = float((img * xs).sum() / total_weight)
    cy = float((img * ys).sum() / total_weight)

    # 重心からのオフセット
    dx = xs - cx
    dy = ys - cy

    # weighted 共分散行列の要素
    cov_xx = float((img * dx * dx).sum() / total_weight)
    cov_xy = float((img * dx * dy).sum() / total_weight)
    cov_yy = float((img * dy * dy).sum() / total_weight)

    cov_matrix = numpy.array([[cov_xx, cov_xy], [cov_xy, cov_yy]], dtype=numpy.float64)

    # 固有値分解（numpy.linalg.eigh は対称行列向けで安定かつ固有値昇順保証）
    eigenvalues, eigenvectors = numpy.linalg.eigh(cov_matrix)
    # eigh は昇順なので，最大固有値は eigenvalues[1]
    primary_eigenvalue = float(eigenvalues[1])
    secondary_eigenvalue = float(eigenvalues[0])

    # 主軸ベクトル（第 1 固有値に対応する固有ベクトル）
    primary_vec = eigenvectors[:, 1]  # shape (2,): [x_component, y_component]
    angle_rad = float(numpy.arctan2(primary_vec[1], primary_vec[0]))
    angle_deg = float(numpy.degrees(angle_rad))
    # -90〜+90 の範囲に折り畳む（主軸は方向性を持たないため）
    if angle_deg > 90.0:
        angle_deg -= 180.0
    elif angle_deg < -90.0:
        angle_deg += 180.0

    # 離心率 eccentricity = √(1 - secondary/primary)
    # primary が 0 の場合（全画素同一値）は 0 を返す
    if primary_eigenvalue <= 0.0:
        eccentricity = 0.0
    else:
        ratio = secondary_eigenvalue / primary_eigenvalue
        # ratio は理論上 [0, 1] だが数値誤差で負になることがあるためクリップ
        eccentricity = float(numpy.sqrt(max(0.0, 1.0 - ratio)))

    return {
        "primary_axis_angle_deg": angle_deg,
        "primary_eigenvalue": primary_eigenvalue,
        "secondary_eigenvalue": secondary_eigenvalue,
        "eccentricity": eccentricity,
    }