pqAutomationApp/app/tests/color_accuracy.py

"""
色准（ΔE2000 / 标准色）相关纯算法。
按 CIE / CIEDE2000 官方定义实现。

标准色度来源：
  - 彩色块(ColorChecker): ColorChecker N Ohta 光谱积分 (CIE1931 2°, D65)，
    已与 Calman 彩色块 Target Y 验证一致（总误差 1.34）。
  - 灰阶: linear（信号值即线性域目标相对亮度），已验证对齐 Calman。
  - 100% 原色/混合色: 由目标色域三原色定义科学推导，按色域分组，
    随信号格式通过 set_active_gamut() 切换。默认 sRGB。

统一规则：
  Target xy = 标准色度
  Target Y  = white_lv × Yr / 100
  Yr 一律是"线性域相对亮度因子"(0~100)，与 xy 同源。
"""

import math
import numpy as np

from app.solid_color_scale import solid_color_rgb as _rgb


D65_X = 0.3127
D65_Y = 0.3290


# ======================================================================
# 彩色块标准 (x, y, Yr)
# 来源：直接抄自 Calman Target xy（X-Rite 厂商标称值），零偏差对齐
# Yr = Calman Target Y / White Target Y(649.5768) × 100
# ======================================================================
_COLOR_PATCH_XYY = {
    "dark skin":     (0.4063, 0.3645, 9.91),
    "light skin":    (0.3780, 0.3562, 35.54),
    "blue sky":      (0.2489, 0.2653, 19.04),
    "foliage":       (0.3416, 0.4319, 13.13),
    "blue flower":   (0.2686, 0.2528, 23.78),
    "bluish green":  (0.2614, 0.3594, 42.42),
    "orange":        (0.5146, 0.4095, 28.59),
    "purplish blue": (0.2147, 0.1891, 11.72),
    "moderate red":  (0.4641, 0.3122, 18.61),
    "purple":        (0.2882, 0.2164, 6.49),
    "yellow green":  (0.3774, 0.4955, 43.66),
    "orange yellow": (0.4749, 0.4427, 42.91),
    "blue":          (0.1883, 0.1349, 6.04),
    "green":         (0.3049, 0.4948, 23.29),
    "red":           (0.5474, 0.3187, 11.57),
    "yellow":        (0.4477, 0.4759, 59.77),
    "magenta":       (0.3738, 0.2440, 18.96),
    "cyan":          (0.2080, 0.2688, 19.71),
}


# ======================================================================
# 100% 原色/混合色：由色域定义科学推导 (x, y, Yr)，按色域分组
# 来源：各色域三原色 + 白点，线性 RGB -> XYZ -> xyY (白场Y=100)
# 注意：sRGB 与 BT.709 共用同一组三原色，数值一致（做别名）
# ======================================================================
_PRIMARY_XYY_BY_GAMUT = {
    "sRGB": {
        "100% Red":     (0.6400, 0.3300, 21.2639),
        "100% Green":   (0.3000, 0.6000, 71.5169),
        "100% Blue":    (0.1500, 0.0600, 7.2192),
        "100% Cyan":    (0.2246, 0.3287, 78.7361),
        "100% Magenta": (0.3209, 0.1542, 28.4831),
        "100% Yellow":  (0.4193, 0.5053, 92.7808),
    },
    "DCI-P3": {
        "100% Red":     (0.6800, 0.3200, 20.9492),
        "100% Green":   (0.2650, 0.6900, 72.1595),
        "100% Blue":    (0.1500, 0.0600, 6.8913),
        "100% Cyan":    (0.2048, 0.3602, 79.0508),
        "100% Magenta": (0.3424, 0.1544, 27.8405),
        "100% Yellow":  (0.4248, 0.5476, 93.1087),
    },
    "BT.2020": {
        "100% Red":     (0.7080, 0.2920, 26.2700),
        "100% Green":   (0.1700, 0.7970, 67.7998),
        "100% Blue":    (0.1310, 0.0460, 5.9302),
        "100% Cyan":    (0.1465, 0.3446, 73.7300),
        "100% Magenta": (0.3682, 0.1471, 32.2002),
        "100% Yellow":  (0.4465, 0.5374, 94.0698),
    },
}
# BT.709 共用 sRGB 三原色，做别名指向，避免重复维护
_PRIMARY_XYY_BY_GAMUT["BT.709"] = _PRIMARY_XYY_BY_GAMUT["sRGB"]

# 当前激活色域（默认 sRGB），随信号格式切换调用 set_active_gamut() 修改
_ACTIVE_GAMUT = "sRGB"


def set_active_gamut(gamut_name):
    """
    切换当前色域（随信号格式选择调用）。
    支持: 'sRGB' / 'BT.709' / 'DCI-P3' / 'BT.2020'。
    """
    global _ACTIVE_GAMUT
    if gamut_name not in _PRIMARY_XYY_BY_GAMUT:
        raise ValueError(
            f"未知色域: {gamut_name}，可选: {list(_PRIMARY_XYY_BY_GAMUT)}"
        )
    _ACTIVE_GAMUT = gamut_name


def get_active_gamut():
    """返回当前激活的色域名。"""
    return _ACTIVE_GAMUT


def _current_primaries():
    """返回当前色域的原色表。"""
    return _PRIMARY_XYY_BY_GAMUT[_ACTIVE_GAMUT]


# ======================================================================
# 测试流程用名 -> 标准来源
#   彩色块: 映射到 _COLOR_PATCH_XYY 的键
#   灰阶  : 映射到 linear 信号值（Yr = 信号值×100）
#   100%原色: 走当前色域的 _PRIMARY_XYY_BY_GAMUT
# ======================================================================
_PATCH_NAME_MAP = {
    "Dark Skin": "dark skin",
    "Light Skin": "light skin",
    "Blue Sky": "blue sky",
    "Foliage": "foliage",
    "Blue Flower": "blue flower",
    "Bluish Green": "bluish green",
    "Orange": "orange",
    "Purplish Blue": "purplish blue",
    "Moderate Red": "moderate red",
    "Purple": "purple",
    "Yellow Green": "yellow green",
    "Orange Yellow": "orange yellow",
    "Blue (Legacy)": "blue",
    "Green (Legacy)": "green",
    "Red (Legacy)": "red",
    "Yellow (Legacy)": "yellow",
    "Magenta (Legacy)": "magenta",
    "Cyan (Legacy)": "cyan",
}

# 灰阶：信号值（线性域目标相对亮度），色度恒为 D65 白点
_GRAYSCALE_SIGNAL = {
    "White": 1.00,
    "Gray 80": 0.80,
    "Gray 65": 0.65,
    "Gray 50": 0.50,
    "Gray 35": 0.35,
    "Black": 0.00,
}


_SDR_COLOR_PATTERNS = [
    ("White", *_rgb(255, 255, 255)),
    ("Gray 80", *_rgb(230, 230, 230)),
    ("Gray 65", *_rgb(209, 209, 209)),
    ("Gray 50", *_rgb(186, 186, 186)),
    ("Gray 35", *_rgb(158, 158, 158)),
    ("Dark Skin", *_rgb(115, 82, 66)),
    ("Light Skin", *_rgb(194, 150, 130)),
    ("Blue Sky", *_rgb(94, 122, 156)),
    ("Foliage", *_rgb(89, 107, 66)),
    ("Blue Flower", *_rgb(130, 128, 176)),
    ("Bluish Green", *_rgb(99, 189, 168)),
    ("Orange", *_rgb(217, 120, 41)),
    ("Purplish Blue", *_rgb(74, 92, 163)),
    ("Moderate Red", *_rgb(194, 84, 97)),
    ("Purple", *_rgb(92, 61, 107)),
    ("Yellow Green", *_rgb(158, 186, 64)),
    ("Orange Yellow", *_rgb(230, 161, 46)),
    ("Blue (Legacy)", *_rgb(51, 61, 150)),
    ("Green (Legacy)", *_rgb(71, 148, 71)),
    ("Red (Legacy)", *_rgb(176, 48, 59)),
    ("Yellow (Legacy)", *_rgb(237, 199, 33)),
    ("Magenta (Legacy)", *_rgb(186, 84, 145)),
    ("Cyan (Legacy)", *_rgb(0, 133, 163)),
    ("100% Red", *_rgb(255, 0, 0)),
    ("100% Green", *_rgb(0, 255, 0)),
    ("100% Blue", *_rgb(0, 0, 255)),
    ("100% Cyan", *_rgb(0, 255, 255)),
    ("100% Magenta", *_rgb(255, 0, 255)),
    ("100% Yellow", *_rgb(255, 255, 0)),
]


# ----------------------------------------------------------------------
# 标准 xy / Yr 解析（统一入口）
# ----------------------------------------------------------------------
def _resolve_reference_xy(name):
    """返回该色块的标准参考 xy。未知则回退 D65。"""
    if name in _GRAYSCALE_SIGNAL:
        return (D65_X, D65_Y)
    if name in _PATCH_NAME_MAP:
        x, y, _ = _COLOR_PATCH_XYY[_PATCH_NAME_MAP[name]]
        return (x, y)
    primaries = _current_primaries()
    if name in primaries:
        x, y, _ = primaries[name]
        return (x, y)
    return (D65_X, D65_Y)


def _resolve_reference_yr(name):
    """
    返回该色块的标准相对亮度因子 Yr (0~100)。
    灰阶: 信号值×100 (linear)。彩色/原色: 光谱或色域 Yr。
    无法确定返回 None。
    """
    if name in _GRAYSCALE_SIGNAL:
        return _GRAYSCALE_SIGNAL[name] * 100.0
    if name in _PATCH_NAME_MAP:
        return _COLOR_PATCH_XYY[_PATCH_NAME_MAP[name]][2]
    primaries = _current_primaries()
    if name in primaries:
        return primaries[name][2]
    return None


def get_target_xyY(name, white_lv):
    """
    返回该色块的完整 target (x, y, Y)。
    Target Y = white_lv × Yr / 100 （统一规则）。
    Yr 不可知时 Y 返回 None，由调用方按实测处理。
    """
    x, y = _resolve_reference_xy(name)
    Yr = _resolve_reference_yr(name)
    if Yr is None or white_lv is None or white_lv <= 0:
        return (x, y, None)
    return (x, y, round(white_lv * Yr / 100.0, 4))


# ----------------------------------------------------------------------
# 兼容旧接口
# ----------------------------------------------------------------------
def get_accuracy_reference_y(name, white_lv):
    """
    返回图表/表格用的参考亮度比例（White=100 缩放）。
    现基于标准 Yr 返回真实比例；无标准 Yr 时回退 100。
    """
    Yr = _resolve_reference_yr(name)
    if Yr is None or white_lv is None or white_lv <= 0:
        return 100.0
    return round(Yr, 4)


# ----------------------------------------------------------------------
# xyY -> XYZ -> Lab（CIE 官方定义）
# ----------------------------------------------------------------------
def _xyY_to_XYZ(x, y, Y):
    """xyY 转 XYZ。y 为 0 时返回全 0，避免除零。"""
    if y <= 0:
        return 0.0, 0.0, 0.0
    X = (x / y) * Y
    Z = ((1.0 - x - y) / y) * Y
    return X, Y, Z


def _XYZ_to_lab(X, Y, Z, Xn, Yn, Zn):
    """
    XYZ 转 CIE L*a*b*。
    Xn, Yn, Zn 为参考白点的绝对 XYZ（Yn 通常为白场亮度 white_lv）。
    """
    delta = 6.0 / 29.0

    def f(t):
        if t > delta ** 3:
            return t ** (1.0 / 3.0)
        return t / (3.0 * delta ** 2) + 4.0 / 29.0

    xr = X / Xn if Xn != 0 else 0.0
    yr = Y / Yn if Yn != 0 else 0.0
    zr = Z / Zn if Zn != 0 else 0.0

    fx, fy, fz = f(xr), f(yr), f(zr)
    L = 116.0 * fy - 16.0
    a = 500.0 * (fx - fy)
    b = 200.0 * (fy - fz)
    return L, a, b


def _xyY_to_lab(x, y, Y, white_x=D65_X, white_y=D65_Y, white_Y=None):
    """
    xyY 直接转 Lab。
    白点默认 D65 色度；white_Y 为白场绝对亮度（用于 L 的归一化基准）。
    若 white_Y 为 None，则退化为以 Y=1 归一化（仅相对比较时可用）。
    """
    if white_Y is None or white_Y <= 0:
        white_Y = 1.0
    X, Y3, Z = _xyY_to_XYZ(x, y, Y)
    Xn, Yn, Zn = _xyY_to_XYZ(white_x, white_y, white_Y)
    return _XYZ_to_lab(X, Y3, Z, Xn, Yn, Zn)


# ----------------------------------------------------------------------
# CIEDE2000（官方公式）
# ----------------------------------------------------------------------
def _delta_e_2000_from_lab(L1, a1, b1, L2, a2, b2, kL=1.0, kC=1.0, kH=1.0):
    C1 = math.sqrt(a1 ** 2 + b1 ** 2)
    C2 = math.sqrt(a2 ** 2 + b2 ** 2)
    C_bar = (C1 + C2) / 2.0

    G = 0.5 * (1 - math.sqrt(C_bar ** 7 / (C_bar ** 7 + 25 ** 7)))

    a1_prime = a1 * (1 + G)
    a2_prime = a2 * (1 + G)

    C1_prime = math.sqrt(a1_prime ** 2 + b1 ** 2)
    C2_prime = math.sqrt(a2_prime ** 2 + b2 ** 2)

    def calc_hue(a_prime, b):
        if a_prime == 0 and b == 0:
            return 0.0
        h = math.degrees(math.atan2(b, a_prime))
        if h < 0:
            h += 360.0
        return h

    h1_prime = calc_hue(a1_prime, b1)
    h2_prime = calc_hue(a2_prime, b2)

    # ΔL', ΔC'
    delta_L_prime = L2 - L1
    delta_C_prime = C2_prime - C1_prime

    # Δh'（官方三分支）
    if C1_prime * C2_prime == 0:
        delta_h_prime = 0.0
    else:
        dh = h2_prime - h1_prime
        if abs(dh) <= 180:
            delta_h_prime = dh
        elif dh > 180:
            delta_h_prime = dh - 360.0
        else:
            delta_h_prime = dh + 360.0

    delta_H_prime = (
        2.0
        * math.sqrt(C1_prime * C2_prime)
        * math.sin(math.radians(delta_h_prime / 2.0))
    )

    # 平均值
    L_bar_prime = (L1 + L2) / 2.0
    C_bar_prime = (C1_prime + C2_prime) / 2.0

    # H_bar'（官方分支）
    if C1_prime * C2_prime == 0:
        H_bar_prime = h1_prime + h2_prime
    else:
        dh_abs = abs(h1_prime - h2_prime)
        if dh_abs <= 180:
            H_bar_prime = (h1_prime + h2_prime) / 2.0
        elif (h1_prime + h2_prime) < 360:
            H_bar_prime = (h1_prime + h2_prime + 360.0) / 2.0
        else:
            H_bar_prime = (h1_prime + h2_prime - 360.0) / 2.0

    T = (
        1
        - 0.17 * math.cos(math.radians(H_bar_prime - 30))
        + 0.24 * math.cos(math.radians(2 * H_bar_prime))
        + 0.32 * math.cos(math.radians(3 * H_bar_prime + 6))
        - 0.20 * math.cos(math.radians(4 * H_bar_prime - 63))
    )

    delta_theta = 30 * math.exp(-(((H_bar_prime - 275) / 25.0) ** 2))
    R_C = 2 * math.sqrt(C_bar_prime ** 7 / (C_bar_prime ** 7 + 25 ** 7))
    R_T = -R_C * math.sin(math.radians(2 * delta_theta))

    S_L = 1 + (0.015 * (L_bar_prime - 50) ** 2) / math.sqrt(
        20 + (L_bar_prime - 50) ** 2
    )
    S_C = 1 + 0.045 * C_bar_prime
    S_H = 1 + 0.015 * C_bar_prime * T

    return math.sqrt(
        (delta_L_prime / (kL * S_L)) ** 2
        + (delta_C_prime / (kC * S_C)) ** 2
        + (delta_H_prime / (kH * S_H)) ** 2
        + R_T * (delta_C_prime / (kC * S_C)) * (delta_H_prime / (kH * S_H))
    )


def calculate_delta_e_2000(
    measured_x,
    measured_y,
    measured_lv,
    standard_x,
    standard_y,
    standard_lv=None,
    white_lv=None,
):
    """
    计算 ΔE 2000 色差。

    Args:
        measured_x, measured_y: 测量的 xy 坐标
        measured_lv: 测量的亮度（cd/m²）
        standard_x, standard_y: 标准的 xy 坐标
        standard_lv: 标准亮度（cd/m²）；默认与 measured_lv 相同
        white_lv: 白场亮度（cd/m²），作为 Lab 的 L 归一化基准；
                  默认取 measured_lv（仅当不传时退化为相对比较）

    Returns:
        float: ΔE 2000 色差值
    """
    if standard_lv is None:
        standard_lv = measured_lv
    if white_lv is None:
        white_lv = measured_lv

    L1, a1, b1 = _xyY_to_lab(measured_x, measured_y, measured_lv, white_Y=white_lv)
    L2, a2, b2 = _xyY_to_lab(standard_x, standard_y, standard_lv, white_Y=white_lv)
    return _delta_e_2000_from_lab(L1, a1, b1, L2, a2, b2)


def calculate_accuracy_delta_e_2000(
    patch_name, measured_x, measured_y, measured_lv, white_lv
):
    """
    色准测试专用 ΔE2000
    标准 xy 来自光谱积分/色域定义（随当前色域）；
    目标 Y 取实测 Y（同亮度下比较色度差异）。
    L 的归一化基准使用白场亮度 white_lv。
    """
    standard_x, standard_y = _resolve_reference_xy(patch_name)
    return calculate_delta_e_2000(
        measured_x,
        measured_y,
        measured_lv,
        standard_x,
        standard_y,
        standard_lv=measured_lv,
        white_lv=white_lv,
    )


def calculate_color_accuracy(measured, standard):
    """计算色差（简化版，xy 欧氏距离 × 1000）"""
    delta_E = {}
    for color in measured.keys():
        dx = measured[color][0] - standard[color][0]
        dy = measured[color][1] - standard[color][1]
        delta_E[color] = np.sqrt(dx * dx + dy * dy) * 1000
    return delta_E


def get_accuracy_color_standards(test_type=None):
    """返回色准标准（patch 名称 -> 参考 xy，随当前色域）。"""
    del test_type
    return {name: _resolve_reference_xy(name) for name, _, _, _ in _SDR_COLOR_PATTERNS}