用Python绘制CIE 1931色度图：从光谱数据到色彩可视化实战
当设计师抱怨显示器色域太窄时，他们其实在谈论什么？为什么专业摄影需要Adobe RGB而不是sRGB？要真正理解这些问题的本质，我们需要回到1931年那个改变色彩科学的坐标系——CIE 1931色度图。本文将带你用Python代码完整复现这个神奇的色彩地图，通过代码理解人眼如何感知颜色边界。
1. 准备工作：理解色彩科学基础
在开始编码之前，我们需要明确几个关键概念。CIE 1931色度图本质上是一个二维投影，它将人眼可见的所有颜色压缩到一个马蹄形区域内。这个马蹄形的边缘被称为光谱轨迹，对应着单一波长的纯色光（从380nm的紫色到780nm的红色）。
注意：色度图的x,y坐标并不直接对应物理波长，而是经过标准观察者函数加权计算的结果
我们需要准备以下Python库：
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import interp1d from matplotlib.patches import Polygon2. 获取标准观察者数据
CIE 1931标准观察者数据是构建色度图的基础，它定义了人眼对可见光谱的响应曲线。我们可以直接从CIE官网获取这些数据，或者使用科学计算库中预置的值：
# CIE 1931 2度标准观察者数据 wavelengths = np.arange(360, 831, 1) # 360-830nm x_bar = np.array([...]) # 填入实际x曲线数据 y_bar = np.array([...]) # 填入实际y曲线数据 z_bar = np.array([...]) # 填入实际z曲线数据为了确保数据准确性，建议使用线性插值处理原始数据：
# 创建插值函数 x_interp = interp1d(wavelengths, x_bar, kind='linear') y_interp = interp1d(wavelengths, y_bar, kind='linear') z_interp = interp1d(wavelengths, z_bar, kind='linear') # 生成更密集的采样点 dense_wl = np.linspace(360, 830, 1000) x_dense = x_interp(dense_wl) y_dense = y_interp(dense_wl) z_dense = z_interp(dense_wl)3. 计算色度坐标与绘制光谱轨迹
色度坐标(x,y)的计算公式看似简单，却蕴含着色彩科学的精髓：
$$ x = \frac{X}{X+Y+Z}, \quad y = \frac{Y}{X+Y+Z} $$
在代码中实现这个转换：
def calculate_xy(X, Y, Z): """计算色度坐标""" sum_xyz = X + Y + Z x = X / sum_xyz y = Y / sum_xyz return x, y # 计算光谱轨迹的色度坐标 spectral_x = [] spectral_y = [] for wl in dense_wl: X = x_interp(wl) Y = y_interp(wl) Z = z_interp(wl) x, y = calculate_xy(X, Y, Z) spectral_x.append(x) spectral_y.append(y)现在可以绘制出著名的马蹄形光谱轨迹：
plt.figure(figsize=(10, 8)) plt.plot(spectral_x, spectral_y, color='black', linewidth=2) plt.title('CIE 1931 Chromaticity Diagram') plt.xlabel('x') plt.ylabel('y') plt.grid(True)4. 添加普朗克轨迹与常见色域
普朗克轨迹（黑体辐射轨迹）展示了不同温度下黑体辐射在色度图上的位置。我们可以通过普朗克公式计算这些点：
def planckian_locus(temperature): """计算指定温度下的黑体辐射色度坐标""" # 实现普朗克公式计算 # 返回x,y坐标 return x, y # 绘制几个典型温度点 temperatures = [2000, 3000, 4000, 5000, 6500, 10000] for temp in temperatures: x, y = planckian_locus(temp) plt.scatter(x, y, color='red') plt.text(x, y, f'{temp}K', fontsize=10)常见色域的三角形边界也是理解色彩空间的关键。以sRGB为例：
# sRGB色域顶点坐标 srgb_red = (0.64, 0.33) srgb_green = (0.30, 0.60) srgb_blue = (0.15, 0.06) # 创建多边形并填充 srgb_triangle = Polygon([srgb_red, srgb_green, srgb_blue], closed=True, fill=True, alpha=0.2, color='blue') plt.gca().add_patch(srgb_triangle)同样方法可以添加Adobe RGB、DCI-P3等其他色域，通过不同颜色和透明度区分它们。
5. 色彩可视化技巧与优化
为了让色度图更直观，我们可以添加以下增强效果：
- 色彩填充：在光谱轨迹内填充近似颜色
from matplotlib.colors import LinearSegmentedColormap # 创建自定义色彩映射 points = np.linspace(0, 1, 100) colors = [...] # 根据色度图定义颜色 cmap = LinearSegmentedColormap.from_list('cie_colormap', colors) # 填充色彩 plt.fill(spectral_x, spectral_y, color='none', edgecolor='none') plt.gca().autoscale_view()- 关键波长标注：标记几个重要波长点
key_wavelengths = [450, 520, 580, 620, 700] # nm for wl in key_wavelengths: idx = np.argmin(np.abs(dense_wl - wl)) plt.scatter(spectral_x[idx], spectral_y[idx], color='black') plt.text(spectral_x[idx], spectral_y[idx], f'{wl}nm', fontsize=10)- 白点标注：标记D65等标准白点
d65 = (0.3127, 0.3290) plt.scatter(d65[0], d65[1], color='white', edgecolor='black', s=100) plt.text(d65[0], d65[1], 'D65', ha='center', va='center')6. 实际应用与扩展思考
完成基础色度图后，我们可以进一步探索几个实用方向：
色域覆盖率计算：比较不同色域的面积
def calculate_area(points): """计算多边形面积""" x = [p[0] for p in points] y = [p[1] for p in points] return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) srgb_area = calculate_area([srgb_red, srgb_green, srgb_blue]) adobe_rgb_area = ... # 类似计算Adobe RGB面积 print(f'sRGB覆盖了{srgb_area/visible_area*100:.1f}%可见色域')色彩转换验证：检查RGB值在色度图中的位置
def rgb_to_xy(r, g, b): """将RGB转换为色度坐标""" # 实现转换矩阵计算 return x, y # 测试纯红色 x, y = rgb_to_xy(1, 0, 0) plt.scatter(x, y, color='red', s=100)色差计算：评估两个颜色的感知差异
def delta_e(x1, y1, x2, y2): """计算两个色度坐标之间的色差""" return np.sqrt((x2-x1)**2 + (y2-y1)**2)在显示器校准项目中，这种可视化特别有用。当我们需要解释为什么某台显示器无法准确再现特定红色时，直接展示该颜色在色度图上的位置与显示器色域边界的关系，比任何口头解释都更有说服力。
