"""放射状輝度プロファイルモジュール．

ROI の中心からの正規化楕円距離を等分割したビンごとに平均輝度を集計し，
中心-周辺方向の輝度変化を連続的に把握するためのプロファイルを生成する．

References:
    [1] MiniTias-Evaluation SPEC_02 「放射状輝度プロファイル」節
        https://github.com/rintoHasegawa/MiniTias-Evaluation
"""

import numpy

from sm_tias_lightsim.evaluation.zones import compute_normalized_distance_map


def calc_radial_profile(
    luminance: numpy.ndarray,
    n_bins: int = 20,
) -> numpy.ndarray:
    """輝度画像から放射状輝度プロファイルを算出する．

    正規化楕円距離 $d \\in [0, 1]$ を `n_bins` 等分のビンに分割し，各ビン内の
    平均輝度を集計する [1]．最外ビンは $d = 1.0$ の画素を欠損させないよう上端を含める．

    Args:
        luminance: 輝度画像（shape: (H, W), dtype: float64）．
        n_bins: ビン数．デフォルトは 20．

    Returns:
        放射状プロファイル配列（shape: (n_bins, 2), dtype: numpy.float64）．
        列 0: ビン中心距離（ビン両端の中点）．
        列 1: そのビンの平均輝度．空ビンは numpy.nan．
    """
    height, width = luminance.shape
    distance_map = compute_normalized_distance_map(height, width)

    # [0, 1] を n_bins 等分するビン境界を生成する
    bin_edges = numpy.linspace(0.0, 1.0, n_bins + 1)
    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0

    # numpy.digitize でビンインデックスを割り当てる（1 始まり, 右端開区間）
    # 最外ビン（d == 1.0）が n_bins + 1 に割り振られないよう right=True を使用する
    bin_indices = numpy.digitize(distance_map.ravel(), bin_edges, right=True)
    flat_luminance = luminance.ravel()

    mean_luminance = numpy.full(n_bins, numpy.nan, dtype=numpy.float64)
    for i in range(1, n_bins + 1):
        mask = bin_indices == i
        if numpy.any(mask):
            mean_luminance[i - 1] = float(numpy.mean(flat_luminance[mask]))
        # 空ビンは初期値 nan のまま維持する

    profile = numpy.column_stack([bin_centers, mean_luminance])
    return profile