别再死记硬背了！用Python+Matplotlib亲手画一遍CIE 1931色度图，理解色彩空间的底层逻辑

用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}")

这个转换带来两个重要特性：

1. 亮度分离 ：Y值单独表示亮度信息
2. 色度纯化 ：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色度图由三个核心几何要素构成：

1. 光谱轨迹 ：单色光在xy平面的投影形成的马蹄形曲线
2. 紫红线 ：连接380nm和780nm的直线，代表非光谱色
3. 白点 ：参考光源(如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")