用Python绘制CIE 1931色度图：从三刺激值到可视化理解的完整实践
当设计师在Pantone色卡前犹豫不决时，当摄影师在Lightroom中调整白平衡时，背后起作用的正是那个神秘的马蹄形图表——CIE 1931色度图。这个看似简单的二维平面，实则是人类视觉感知的数学映射，它定义了"物理光信号"如何转化为"心理颜色感知"的桥梁。本文将带您用Python代码亲手构建这个色彩科学的基石，在数据计算与图形绘制中理解色彩空间的本质。
1. 色彩科学基础与实验数据准备
CIE 1931系统的核心在于将无限可能的可见光谱压缩为三个数值——X、Y、Z三刺激值。这源于1920年代Wright和Guild的突破性发现：人类视网膜上的三种视锥细胞对光的响应曲线，可以通过特定比例的红、绿、蓝光刺激来模拟。
获取标准观察者数据 是第一步。国际照明委员会(CIE)提供的标准数据包含两个关键部分：
import numpy as np
# CIE 1931 2°标准观察者数据 (波长间隔5nm)
wavelengths = np.arange(380, 781, 5) # 380-780nm可见光谱
cmf_1931 = np.array([
[0.0014, 0.0000, 0.0065], # 380nm
[0.0022, 0.0001, 0.0105], # 385nm
... # 完整数据应包含81组值
[0.0003, 0.0000, 0.0001] # 780nm
])
注意：实际应用中建议使用CIE发布的完整数据文件，上述仅为示例格式。X、Y、Z分量分别对应人眼对长波、中波、短波的敏感度曲线。
理解这些数据的物理意义至关重要：
X分量 ：模拟L型视锥细胞对红光的响应
Y分量 ：同时代表M型视锥细胞响应和亮度感知
Z分量 ：对应S型视锥细胞的蓝光敏感度
2. 从光谱到色度坐标的计算路径
色度坐标的计算本质上是将三维XYZ空间投影到二维平面。这个降维过程通过归一化实现：
def xyz_to_xy(X, Y, Z):
"""将XYZ三刺激值转换为xy色度坐标"""
sum_xyz = X + Y + Z
return X/sum_xyz, Y/sum_xyz
# 示例：计算D65标准光源的色度坐标
D65_XYZ = [95.047, 100.00, 108.883] # 标准D65光源
x_D65, y_D65 = xyz_to_xy(*D65_XYZ)
print(f"D65色度坐标: x={x_D65:.4f}, y={y_D65:.4f}")
这个转换带来两个重要特性：
亮度分离 ：Y值单独表示亮度信息
色度纯化 ：xy坐标仅反映颜色的"色调"和"饱和度"
光谱轨迹计算 是绘制色度图的关键步骤。我们需要对可见光谱的每个波长点进行计算：
spectral_xy = []
for i in range(len(wavelengths)):
X, Y, Z = cmf_1931[i]
x, y = xyz_to_xy(X, Y, Z)
spectral_xy.append((x, y))
3. 构建色度图的几何要素
完整的CIE 1931色度图由三个核心几何要素构成：
光谱轨迹 ：单色光在xy平面的投影形成的马蹄形曲线
紫红线 ：连接380nm和780nm的直线，代表非光谱色
白点 ：参考光源(如D65)的位置
Matplotlib绘制实现 ：
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
fig, ax = plt.subplots(figsize=(10, 8))
# 绘制光谱轨迹
spectral_x, spectral_y = zip(*spectral_xy)
ax.plot(spectral_x, spectral_y, color='black', linewidth=2)
# 绘制紫红线
ax.plot([spectral_x[0], spectral_x[-1]],
[spectral_y[0], spectral_y[-1]],
color='purple', linestyle='--')
# 标记白点
ax.scatter(x_D65, y_D65, color='white', edgecolor='black', s=100, label='D65')
# 填充色域
spectrum_polygon = Polygon(list(zip(spectral_x, spectral_y)) +
[(spectral_x[-1], spectral_y[-1]),
(spectral_x[0], spectral_y[0])],
closed=True, alpha=0.3)
ax.add_patch(spectrum_polygon)
4. 色度图的高级可视化技巧
基础图形绘制完成后，我们可以通过以下增强手段提升信息传达效果：
色彩填充技术 ：
from matplotlib.colors import LinearSegmentedColormap
# 创建色度图专用的非均匀色彩映射
cie_cmap = LinearSegmentedColormap.from_list('cie_map',
[(0, (1,0,0)), (0.3, (1,1,0)), (0.6, (0,1,0)),
(0.8, (0,1,1)), (1, (0,0,1))])
# 在色域多边形内创建网格点进行颜色插值
xx, yy = np.meshgrid(np.linspace(0, 0.8, 200),
np.linspace(0, 0.9, 200))
inside = []
for x, y in zip(xx.flatten(), yy.flatten()):
if spectrum_polygon.contains_point((x, y)):
inside.append((x, y))
等色调线与等饱和度线 ：
# 计算从白点到光谱轨迹的等色调线
for hue_angle in np.linspace(0, 360, 12):
rad = np.deg2rad(hue_angle)
end_x = x_D65 + 0.5 * np.cos(rad)
end_y = y_D65 + 0.5 * np.sin(rad)
ax.plot([x_D65, end_x], [y_D65, end_y],
color='gray', alpha=0.3, linestyle=':')
关键标注元素 ：
# 标记主要波长点
key_wavelengths = [450, 500, 550, 580, 600, 650]
for wl in key_wavelengths:
idx = np.where(wavelengths == wl)[0][0]
ax.annotate(f'{wl}nm', (spectral_x[idx], spectral_y[idx]),
textcoords="offset points", xytext=(5,5),
ha='left', fontsize=10)
5. 色度图的实际应用解析
理解色度图的最终目的是应用。以下是三个典型场景：
1. 色域比较 ：
# 定义sRGB和Adobe RGB色域边界
srgb_vertices = [(0.64,0.33), (0.30,0.60), (0.15,0.06)]
adobe_rgb_vertices = [(0.64,0.33), (0.21,0.71), (0.15,0.06)]
# 绘制色域对比
ax.add_patch(Polygon(srgb_vertices, closed=True,
fill=False, edgecolor='blue', linewidth=2))
ax.add_patch(Polygon(adobe_rgb_vertices, closed=True,
fill=False, edgecolor='green', linewidth=2))
2. 颜色混合预测 ：
def mix_colors(xy1, xy2, ratio):
"""预测两种颜色混合后的色度坐标"""
mix_x = ratio*xy1[0] + (1-ratio)*xy2[0]
mix_y = ratio*xy1[1] + (1-ratio)*xy2[1]
return (mix_x, mix_y)
red = (0.64, 0.33)
green = (0.30, 0.60)
for ratio in np.linspace(0, 1, 10):
mix = mix_colors(red, green, ratio)
ax.scatter(*mix, color='black', s=20)
3. 色温轨迹可视化 ：
# 普朗克轨迹数据 (黑体辐射色温线)
planckian_locus = [(0.4476,0.4074), (0.3484,0.3516),
(0.3101,0.3162), (0.2856,0.2827)]
ax.plot(*zip(*planckian_locus), color='red',
linestyle='-.', label='Planckian Locus')
6. 工程实践中的注意事项
在实际项目中应用色度图时，有几个关键点需要特别注意：
色彩转换的精度问题 ：
8位RGB到XYZ的转换会引入量化误差
不同白点假设会导致色度坐标偏移
显示器色域限制可能使理论计算无法实现
def rgb_to_xyz(rgb, gamma=2.2):
"""考虑gamma校正的RGB转XYZ"""
linear_rgb = np.where(rgb <= 0.04045,
rgb/12.92,
((rgb+0.055)/1.055)**gamma)
# 此处应有完整的转换矩阵计算
return xyz
视觉均匀性问题 ： CIE 1931的一个主要缺陷是色度图上的距离与视觉感知差异不一致。这导致了后续CIELUV等均匀色彩空间的开发：
def xy_to_Luv(x, y, Y=100):
"""将xyY转换到CIELUV空间"""
u_prime = 4*x / (-2*x + 12*y + 3)
v_prime = 9*y / (-2*x + 12*y + 3)
# 完整计算还应考虑参考白点等参数
return (116*(Y/100)**(1/3)-16, 13*L*(u_prime-u_ref), 13*L*(v_prime-v_ref))
现代色彩管理实践 ：
使用ICC配置文件确保跨设备一致性
在图像处理管线中保持线性色彩空间计算
对广色域内容进行适当的元数据标记
import PIL.ImageCms
# 创建从sRGB到Adobe RGB的色彩转换管道
srgb_profile = PIL.ImageCms.createProfile("sRGB")
adobe_rgb_profile = PIL.ImageCms.createProfile("Adobe RGB")
transform = PIL.ImageCms.buildTransform(
srgb_profile, adobe_rgb_profile, "RGB", "RGB")