从色卡到代码：手把手教你用Python实现CIE 1931色度图转换（附完整代码）

从色卡到代码：手把手教你用Python实现CIE 1931色度图转换

在数字图像处理领域，色彩科学是连接物理世界与数字表达的桥梁。当我们谈论"苹果红"或"天空蓝"时，如何将这些主观感受转化为计算机可处理的客观数据？这正是CIE 1931色度图要解决的核心问题。作为国际色彩标准的基石，这套系统通过数学建模将人眼对颜色的感知量化，为现代显示器校准、印刷色彩管理、影视调色等应用提供了统一语言。

本文将带您深入色彩科学的实践层面，从基本原理到完整代码实现。不同于教科书式的理论阐述，我们会通过Python代码逐步构建色度图转换工具，让抽象概念变得可触摸。您将学到：

1. 如何用NumPy处理光谱数据
2. 颜色匹配函数的实际应用
3. XYZ坐标系的工程意义
4. 可视化色度图的技巧

1. 色彩科学基础与环境准备

1.1 理解CIE 1931的核心概念

CIE 1931系统的精妙之处在于它用三个虚构的原色X、Y、Z建立了无负值的色彩空间。这种设计解决了实际物理原色（如RGB）在匹配某些光谱色时出现负值的尴尬。三个关键特性使其成为工业标准：

• Y通道代表亮度：直接对应人眼对绿光的敏感曲线
• 色度坐标归一化：x = X/(X+Y+Z), y = Y/(X+Y+Z)
• 马蹄形色域边界：包含人眼可见的所有颜色

实验数据显示，人眼对约555nm的绿光最敏感，这解释了为什么Y通道被赋予亮度信息。下表对比了不同色彩空间的特点：

1.2 配置Python环境

推荐使用Anaconda创建专用环境：

conda create -n color_science python=3.9 conda activate color_science pip install numpy matplotlib scipy pandas

关键库的作用：

• NumPy：处理光谱数据矩阵运算
• Matplotlib：绘制色度图及可视化
• SciPy：插值计算等科学计算
• Pandas：整理CIE标准观察者数据

提示：建议使用Jupyter Notebook进行交互式开发，方便实时查看色度图生成效果

2. 加载与处理CIE标准数据

2.1 获取标准观察者数据

CIE官方提供了标准观察者的光谱三刺激值，这是整个系统的基石。我们可以直接从CSV文件加载：

import pandas as pd # 加载CIE 1931 2°标准观察者数据 df = pd.read_csv('cie_1931_2deg.csv', header=None, names=['wavelength', 'x_bar', 'y_bar', 'z_bar']) # 提取380nm-780nm范围（5nm间隔） wavelengths = df['wavelength'].values x_bar = df['x_bar'].values y_bar = df['y_bar'].values z_bar = df['z_bar'].values

数据格式说明：

• 每行对应特定波长的三刺激值
• 波长范围通常为380-780nm（可见光谱）
• x̄, ȳ, z̄代表单位能量的光谱色匹配函数

2.2 数据预处理技巧

实际应用中常需要处理不完整或非常规间隔的数据：

from scipy.interpolate import interp1d # 创建插值函数 interp_x = interp1d(wavelengths, x_bar, kind='cubic') interp_y = interp1d(wavelengths, y_bar, kind='cubic') interp_z = interp1d(wavelengths, z_bar, kind='cubic') # 生成1nm间隔的密集数据 dense_wl = np.arange(380, 781) dense_x = interp_x(dense_wl) dense_y = interp_y(dense_wl) dense_z = interp_z(dense_wl)

注意：CIE数据通常以5nm或10nm为间隔，插值可获得更平滑的色度图边界

3. 实现光谱到XYZ的转换

3.1 光谱功率分布的处理

假设我们有一个LED的光谱功率分布(SPD)测量数据：

# 示例：白色LED的SPD (任意单位) led_spd = np.array([ [400, 0.02], [450, 0.15], [500, 0.35], [550, 0.8], [600, 0.7], [650, 0.4], [700, 0.1] ]) # 归一化处理 led_wl = led_spd[:, 0] led_power = led_spd[:, 1] / np.max(led_spd[:, 1])

3.2 计算三刺激值

根据CIE定义，XYZ三刺激值是SPD与匹配函数的乘积积分：

def spectrum_to_xyz(wavelengths, power, x, y, z, ref_wl): """将光谱转换为XYZ值""" # 插值使SPD与标准观察者数据对齐 interp_power = interp1d(wavelengths, power, bounds_error=False, fill_value=0) aligned_power = interp_power(ref_wl) # 数值积分（黎曼和） delta = ref_wl[1] - ref_wl[0] X = np.sum(aligned_power * x) * delta Y = np.sum(aligned_power * y) * delta Z = np.sum(aligned_power * z) * delta return X, Y, Z # 计算示例LED的XYZ led_XYZ = spectrum_to_xyz(led_wl, led_power, dense_x, dense_y, dense_z, dense_wl) print(f"LED XYZ值: {led_XYZ}")

3.3 色度坐标计算

将XYZ转换为色度坐标(x,y)：

def xyz_to_xy(X, Y, Z): """XYZ转色度坐标""" sum_XYZ = X + Y + Z x = X / sum_XYZ y = Y / sum_XYZ return x, y led_x, led_y = xyz_to_xy(*led_XYZ) print(f"色度坐标: x={led_x:.4f}, y={led_y:.4f}")

4. 绘制CIE 1931色度图

4.1 构建色度图边界

色度图的马蹄形边界由单色光谱的色度坐标构成：

def calculate_spectrum_locus(x, y, z): """计算光谱轨迹""" sum_xyz = x + y + z x_chroma = x / sum_xyz y_chroma = y / sum_xyz return x_chroma, y_chroma # 计算光谱轨迹 spectrum_x, spectrum_y = calculate_spectrum_locus(dense_x, dense_y, dense_z) # 添加紫红线（连接380nm和780nm） purple_x = np.linspace(spectrum_x[0], spectrum_x[-1], 50) purple_y = np.linspace(spectrum_y[0], spectrum_y[-1], 50)

4.2 可视化实现

使用Matplotlib创建专业级色度图：

import matplotlib.pyplot as plt from matplotlib.patches import Polygon plt.figure(figsize=(10, 8)) # 创建色度图多边形 xy_points = np.column_stack([np.concatenate([spectrum_x, purple_x[::-1]]), np.concatenate([spectrum_y, purple_y[::-1]])]) poly = Polygon(xy_points, closed=True, alpha=0.3) plt.gca().add_patch(poly) # 绘制光谱轨迹 plt.plot(spectrum_x, spectrum_y, color='black', linewidth=1.5) # 标记示例点 plt.scatter([led_x], [led_y], color='red', s=100, label='LED光源') plt.text(led_x, led_y, f' ({led_x:.3f}, {led_y:.3f})', fontsize=10) # 添加标准白点 (D65) plt.scatter([0.3127], [0.3290], color='white', edgecolor='black', s=100, label='D65') plt.text(0.3127, 0.3290, ' D65', fontsize=10) plt.xlabel('x chromaticity') plt.ylabel('y chromaticity') plt.title('CIE 1931 Chromaticity Diagram') plt.grid(True, linestyle='--', alpha=0.5) plt.legend() plt.xlim(0, 0.8) plt.ylim(0, 0.9) plt.show()

4.3 高级功能扩展

为色度图添加实用功能：

# 色域覆盖计算 def calculate_gamut_coverage(points_x, points_y, ref_x, ref_y): """计算自定义色域在CIE色度图中的覆盖面积""" from scipy.spatial import ConvexHull target_hull = ConvexHull(np.column_stack([points_x, points_y])) ref_hull = ConvexHull(np.column_stack([ref_x, ref_y])) return target_hull.volume / ref_hull.volume # 示例：计算sRGB色域覆盖率 srgb_x = [0.64, 0.30, 0.15] srgb_y = [0.33, 0.60, 0.06] coverage = calculate_gamut_coverage(srgb_x, srgb_y, spectrum_x, spectrum_y) print(f"sRGB色域覆盖率: {coverage*100:.1f}%")

5. 实际应用案例分析

5.1 显示器色域验证

检查设计稿颜色是否超出显示器色域：

def check_color_in_gamut(x, y, gamut_x, gamut_y): """检查颜色点是否在目标色域内""" from matplotlib.path import Path gamut_path = Path(np.column_stack([gamut_x, gamut_y])) return gamut_path.contains_point((x, y)) # 示例：检查P3色域 p3_x = [0.680, 0.265, 0.150] p3_y = [0.320, 0.690, 0.060] design_color = (0.5, 0.3) # 设计稿某颜色 if check_color_in_gamut(*design_color, p3_x, p3_y): print("颜色在P3色域内") else: print("警告：颜色超出P3色域")

5.2 色彩转换工具开发

构建RGB到XYZ的转换管道：

def srgb_to_xyz(r, g, b): """sRGB转CIE XYZ""" # 线性化 def inverse_gamma(u): return u/12.92 if u <= 0.04045 else ((u+0.055)/1.055)**2.4 r_lin = inverse_gamma(r) g_lin = inverse_gamma(g) b_lin = inverse_gamma(b) # 转换矩阵 (sRGB转XYZ) M = np.array([ [0.4124564, 0.3575761, 0.1804375], [0.2126729, 0.7151522, 0.0721750], [0.0193339, 0.1191920, 0.9503041] ]) xyz = M @ np.array([r_lin, g_lin, b_lin]) return xyz # 示例：转换苹果红 apple_red = (0.8, 0.1, 0.1) # sRGB值 apple_xyz = srgb_to_xyz(*apple_red) apple_xy = xyz_to_xy(*apple_xyz) print(f"苹果红的色度坐标: {apple_xy}")

6. 性能优化与工程实践

6.1 加速计算的技巧

处理大批量颜色时，向量化操作可显著提升性能：

def batch_spectrum_to_xyz(spectra, ref_wl, x, y, z): """批量处理光谱数据""" delta = ref_wl[1] - ref_wl[0] xyz = np.zeros((len(spectra), 3)) for i, (wl, power) in enumerate(spectra): interp_p = interp1d(wl, power, bounds_error=False, fill_value=0)(ref_wl) xyz[i, 0] = np.sum(interp_p * x) * delta xyz[i, 1] = np.sum(interp_p * y) * delta xyz[i, 2] = np.sum(interp_p * z) * delta return xyz # 使用示例 multiple_leds = [ (np.array([450, 550, 650]), np.array([0.2, 0.7, 0.3])), # LED A (np.array([470, 520, 620]), np.array([0.1, 0.8, 0.2])) # LED B ] batch_xyz = batch_spectrum_to_xyz(multiple_leds, dense_wl, dense_x, dense_y, dense_z)

6.2 常见问题排查

调试色彩转换时的实用检查点：

1. 数据范围验证：

   • 确认输入光谱波长在380-780nm范围内
   • 检查XYZ值是否均为正数

2. 归一化问题：

   # 验证色度坐标总和 sum_xy = x + y + (1 - x - y) # 应≈1

3. 可视化诊断：

   # 绘制SPD与匹配函数的叠加图 plt.plot(dense_wl, dense_y, label='ȳ') plt.plot(led_wl, led_power, label='LED SPD') plt.legend()

在开发色彩管理工具时，曾遇到因插值方法不当导致色度坐标偏移的问题。改用三次样条插值后，色度图边界平滑度显著提升，特别是在490-500nm的青色区域。另一个经验是始终验证白点的位置——D65应准确落在(0.3127, 0.3290)，这是判断整个转换管道是否正确的重要基准。