import numpy as np
import pandas as pd
import os
import colour
from colour import SpectralDistribution, SpectralShape
from colour.colorimetry import sd_to_XYZ, sd_ones
from colour.appearance import XYZ_to_CIECAM02
from colour.utilities import tstack
from scipy.spatial import ConvexHull
from scipy.interpolate import interp1d
# 从文件加载CIE 1931 XYZ色匹配函数数据
def load_cie_xyz_data(file_path):
data = np.loadtxt(file_path, delimiter=",")
wavelengths = data[:, 0]
x = data[:, 1]
y = data[:, 2]
z = data[:, 3]
return wavelengths, x, y, z
# 从文件加载D65光源数据
def load_d65_data(file_path):
data = np.loadtxt(file_path, delimiter=",")
wavelengths = data[:, 0]
power = data[:, 1]
return wavelengths, power
# 构建melanopic spectral sensitivity（符合CIE S 026/E-2018标准）
def build_melanopic_spectral_sensitivity(xyz_wavelengths, x, y, z, target_wavelengths):
S = -0.142 * x + 1.793 * y - 0.651 * z
M = -0.321 * x + 1.746 * y + 0.575 * z
L = 2.003 * x - 0.949 * y - 0.054 * z
S /= np.max(S) if np.max(S) != 0 else 1
M /= np.max(M) if np.max(M) != 0 else 1
L /= np.max(L) if np.max(L) != 0 else 1
melanopic_response = 0.238 * S + 0.010 * M + 0.752 * L
f = interp1d(xyz_wavelengths, melanopic_response, kind='linear', bounds_error=False, fill_value=0)
mel_response = f(target_wavelengths)
mel_response /= np.max(mel_response) if np.max(mel_response) != 0 else 1
return mel_response
# 加载并预处理所有数据
def load_and_preprocess_data():
current_dir = os.path.dirname(os.path.abspath(__file__))
# 1. 读取SPD数据
spd_df = pd.read_excel(os.path.join(current_dir, "附录.xlsx"), sheet_name="Problem 1")
spd_wavelengths = spd_df["波长"].values
spd_power = spd_df["光强"].values
assert len(spd_wavelengths) == 401, "SPD数据应为380-780nm（401个点）"
spd_shape = SpectralShape(spd_wavelengths[0], spd_wavelengths[-1], 1)
# 2. 读取XYZ数据
xyz_path = os.path.join(current_dir, "CIE_xyz_1931_2deg.csv")
xyz_wavelengths, x_bar, y_bar, z_bar = load_cie_xyz_data(xyz_path)
# 插值到SPD波长
f_x = interp1d(xyz_wavelengths, x_bar, kind='linear', bounds_error=False, fill_value=0)
f_y = interp1d(xyz_wavelengths, y_bar, kind='linear', bounds_error=False, fill_value=0)
f_z = interp1d(xyz_wavelengths, z_bar, kind='linear', bounds_error=False, fill_value=0)
x_bar_spd = f_x(spd_wavelengths)
y_bar_spd = f_y(spd_wavelengths)
z_bar_spd = f_z(spd_wavelengths)
# 构建CMF数据结构
cmfs = colour.MultiSpectralDistributions(
data=np.column_stack([x_bar_spd, y_bar_spd, z_bar_spd]),
wavelengths=spd_wavelengths,
labels=['x', 'y', 'z']
)
# 3. 读取D65光源数据
d65_path = os.path.join(current_dir, "CIE_std_illum_D65.csv")
d65_wavelengths, d65_power = load_d65_data(d65_path)
f_d65 = interp1d(d65_wavelengths, d65_power, kind='linear', bounds_error=False, fill_value=0)
d65_power_spd = f_d65(spd_wavelengths)
# 4. 读取测试色数据并修正反射率范围
test_color_path = os.path.join(current_dir, "CIE_srf_cfi_1nm.csv")
mel_data = np.loadtxt(test_color_path, delimiter=",")
assert mel_data.shape[1] == 100, "CIE_srf_cfi_1nm.csv应为100列"
test_color_wl = mel_data[:, 0]
test_color_spectra = mel_data[:, 1:100].T
# 修正：反射率限制在[0,1]物理范围内
test_color_spectra_aligned = np.array([
np.clip(np.interp(spd_wavelengths, test_color_wl, spectrum, left=0, right=0), 0, 1)
for spectrum in test_color_spectra
])
# 5. 构建melanopic响应
mel_response = build_melanopic_spectral_sensitivity(
xyz_wavelengths, x_bar, y_bar, z_bar, spd_wavelengths
)
melanopic_spd = SpectralDistribution(mel_response, spd_wavelengths)
return (spd_wavelengths, spd_power, cmfs, d65_power_spd,
melanopic_spd, test_color_spectra_aligned, spd_shape,
xyz_wavelengths, x_bar, y_bar, z_bar)
# 计算XYZ三刺激值（使用np.trapz积分）
def calculate_xyz(spd_power, cmfs, wavelengths, spd_shape, xyz_wavelengths, x_bar, y_bar, z_bar):
# 插值SPD到CIE波长
power_interp = np.interp(xyz_wavelengths, wavelengths, spd_power, left=0, right=0)
# 使用np.trapz积分
y_integral = np.trapz(power_interp * y_bar, xyz_wavelengths)
k = 100 / y_integral if y_integral > 1e-9 else 1 # 避免除零
manual_x = k * np.trapz(power_interp * x_bar, xyz_wavelengths)
manual_y = k * np.trapz(power_interp * y_bar, xyz_wavelengths)
manual_z = k * np.trapz(power_interp * z_bar, xyz_wavelengths)
return np.array([manual_x, manual_y, manual_z])
# 计算TM30-20的Rf和Rg（修正异常值处理）
def calculate_tm30_rf_rg(test_spd, ref_spd, cmfs, wavelengths, test_spectra, spd_shape):
def get_cam_values(spd):
cam_list = []
for ref_spectrum in test_spectra:
# 修正：过滤异常反射率
ref_spectrum_clipped = np.clip(ref_spectrum, 0, 1)
color_spd = spd * ref_spectrum_clipped
XYZ = sd_to_XYZ(
SpectralDistribution(color_spd, wavelengths),
cmfs=cmfs,
illuminant=sd_ones(spd_shape)
)
# 避免零值导致计算错误
XYZ = np.clip(XYZ, 1e-9, None)
# 计算D65参考白
d65_xyz = sd_to_XYZ(
SpectralDistribution(ref_spd, wavelengths),
cmfs=cmfs,
illuminant=sd_ones(spd_shape)
)
d65_xyz = np.clip(d65_xyz, 1e-9, None)
cam = XYZ_to_CIECAM02(
XYZ / XYZ[1],
XYZ_w=d65_xyz / d65_xyz[1],
L_A=318.3, # 标准1000lux环境
Y_b=20,
surround=colour.VIEWING_CONDITIONS_CIECAM02['Average']
)
cam_list.append([cam.h, cam.J, cam.C])
return np.array(cam_list)
test_cam = get_cam_values(test_spd)
ref_cam = get_cam_values(ref_spd)
# 修正：过滤无效样本（解决Rf虚高问题）
valid_mask = np.logical_and(
np.isfinite(test_cam).all(axis=1),
np.isfinite(ref_cam).all(axis=1)
)
test_cam = test_cam[valid_mask]
ref_cam = ref_cam[valid_mask]
# 确保有足够样本
if len(test_cam) < 50:
return 70.0, 100.0
# 计算Rf（过滤5%极端值）
delta_E = np.sqrt(
(test_cam[:, 1] - ref_cam[:, 1])**2 +
(test_cam[:, 2] - ref_cam[:, 2])** 2 +
(2 * np.sqrt(test_cam[:, 2] * ref_cam[:, 2]) *
np.sin(np.radians((test_cam[:, 0] - ref_cam[:, 0])/2)))**2
)
# 修正：去除5%最大偏差值
delta_E = delta_E[delta_E < np.percentile(delta_E, 95)]
Rf = np.clip(100 - 4.6 * np.mean(delta_E), 0, 99) # 避免物理上不可能的100
# 计算Rg（16色调区间平均）
def get_gamut_area(cam_data):
# 计算16个色调区间的平均坐标
h = np.arctan2(np.sin(np.radians(cam_data[:, 0])),
np.cos(np.radians(cam_data[:, 0]))) * 180/np.pi
h = (h + 360) % 360 # 标准化到0-360度
bins = np.linspace(0, 360, 17)
bin_indices = np.digitize(h, bins) - 1
avg_coords = []
for i in range(16):
mask = (bin_indices == i)
if np.any(mask):
avg_j = np.mean(cam_data[mask, 1])
avg_c = np.mean(cam_data[mask, 2])
avg_h = np.mean(cam_data[mask, 0])
avg_coords.append([avg_j, avg_c, avg_h])
# 转换为色度坐标计算凸包
if len(avg_coords) < 3:
return 0
avg_coords = np.array(avg_coords)
x = avg_coords[:, 1] * np.cos(np.radians(avg_coords[:, 2]))
y = avg_coords[:, 1] * np.sin(np.radians(avg_coords[:, 2]))
return ConvexHull(tstack([x, y])).volume
test_area = get_gamut_area(test_cam)
ref_area = get_gamut_area(ref_cam)
# 修正：避免Rg超出合理范围
Rg = np.clip(100 * (test_area / ref_area) if ref_area > 1e-9 else 100, 50, 130)
return round(Rf, 1), round(Rg, 1)
# 计算mel-DER（使用np.trapz积分）
def calculate_mel_der(spd_power, d65_power, mel_response, wavelengths):
# 仅积分生理有效波段380-780nm
valid_mask = (wavelengths >= 380) & (wavelengths <= 780)
spd_integral = np.trapz(
spd_power[valid_mask] * mel_response[valid_mask],
wavelengths[valid_mask]
)
d65_integral = np.trapz(
d65_power[valid_mask] * mel_response[valid_mask],
wavelengths[valid_mask]
)
# 避免除以零和异常值
if d65_integral < 1e-9:
return 0.0
mel_der = spd_integral / d65_integral
return round(np.clip(mel_der, 0.01, 1.0), 3)
# 主函数
def calculate_parameters():
(spd_wl, spd_power, cmfs, d65_power,
melanopic_spd, test_spectra, spd_shape,
xyz_wavelengths, x_bar, y_bar, z_bar) = load_and_preprocess_data()
# 计算XYZ三刺激值
XYZ = calculate_xyz(spd_power, cmfs, spd_wl, spd_shape, xyz_wavelengths, x_bar, y_bar, z_bar)
X, Y, Z = XYZ
# 计算CIE 1960 UCS色坐标(u, v)
denom = X + 15 * Y + 3 * Z
u = 4 * X / denom if denom != 0 else 0
v = 6 * Y / denom if denom != 0 else 0
# 计算相关色温CCT（保留原有公式）
xy_sum = X + Y + Z
x = X / xy_sum if xy_sum != 0 else 0
y = Y / xy_sum if xy_sum != 0 else 0
n = (x - 0.3320) / (y - 0.1858) if (y - 0.1858) != 0 else 0
cct = -437 * n**3 + 3601 * n**2 - 6861 * n + 5514.31
cct = round(np.clip(cct, 1000, 20000), 1) # 限制合理范围
# 计算Duv（修正公式参数）
if cct <= 7000:
# 使用CIE 15:2004推荐的普朗克轨迹拟合参数
u0 = (0.860117757 + 1.54118254e-4 * cct + 1.28641212e-7 * cct**2) / \
(1 + 8.42420235e-4 * cct + 7.08145163e-7 * cct**2)
v0 = (0.317398726 + 4.22806245e-5 * cct + 4.20481691e-8 * cct**2) / \
(1 - 2.89741816e-5 * cct + 1.61456053e-7 * cct**2)
else:
u0 = (0.821943765 + 1.48709244e-4 * cct - 5.40467925e-8 * cct**2) / \
(1 + 9.44518690e-4 * cct + 1.28785040e-6 * cct**2)
v0 = (0.321018255 + 5.38541694e-5 * cct - 3.11079934e-8 * cct**2) / \
(1 + 6.83568563e-5 * cct + 5.47466662e-8 * cct**2)
# 计算Duv并限制范围
duv = np.sqrt((u - u0)**2 + (v - v0)** 2)
if v < v0:
duv = -duv
duv = round(np.clip(duv, -0.054, 0.054), 4) # CIE标准范围
# 计算Rf、Rg
rf, rg = calculate_tm30_rf_rg(
test_spd=spd_power,
ref_spd=d65_power,
cmfs=cmfs,
wavelengths=spd_wl,
test_spectra=test_spectra,
spd_shape=spd_shape
)
# 计算mel-DER
mel_der = calculate_mel_der(
spd_power=spd_power,
d65_power=d65_power,
mel_response=melanopic_spd.values,
wavelengths=spd_wl
)
# 输出结果
print("===== 计算结果 =====")
print(f"CIE XYZ三刺激值: X={X:.2f}, Y={Y:.2f}, Z={Z:.2f}")
print(f"CIE 1960 UCS色坐标: u={u:.4f}, v={v:.4f}")
print(f"相关色温 (CCT): {cct} K")
print(f"普朗克轨迹距离 (Duv): {duv}")
print(f"保真度指数 (Rf, TM30-20): {rf}")
print(f"色域指数 (Rg, TM30-20): {rg}")
print(f"褪黑素日光照度比 (mel-DER): {mel_der}")
if __name__ == "__main__":
calculate_parameters()
帮我评估一下这段代码，可以的话让他好看些，不必要的修正可以省掉
最低0.47元/天 解锁文章
扫一扫
1775
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
