别再只调RGB了！用Python+Matplotlib可视化CIE色度图，理解色彩管理的底层逻辑

用Python绘制CIE色度图：从色彩理论到代码实践

色彩科学是现代数字视觉领域的基石，而CIE色度图则是理解这一复杂系统的关键工具。作为开发者，我们常常需要处理色彩空间转换、色域映射和设备校准等问题，但很少有人真正深入探究背后的数学原理。本文将带你用Python和Matplotlib从零构建CIE色度图，通过代码实现揭示色彩管理的底层逻辑。

1. CIE色度系统基础与Python实现

CIE 1931 XYZ色彩空间是人类视觉研究的里程碑式成果。它的核心思想是将所有可见光颜色表示为三个虚构的刺激值X、Y、Z。有趣的是，这些值并非随意设定——Y分量直接对应人眼感知的亮度，而X和Z则代表色度信息。

在Python中，我们可以用numpy数组来表示这些值：

import numpy as np

# 标准观察者数据 - CIE 1931 2度视场
lambda_nm = np.arange(360, 831, 1)  # 波长范围360-830nm
x_bar = np.loadtxt('cie_xyz_1931_2deg_x.txt')  # X三刺激函数
y_bar = np.loadtxt('cie_xyz_1931_2deg_y.txt')  # Y三刺激函数(亮度)
z_bar = np.loadtxt('cie_xyz_1931_2deg_z.txt')  # Z三刺激函数

注意：实际应用中需要加载标准CIE观察者数据文件，这些数据可以从CIE官网或色彩科学库中获取

色度坐标的计算公式看似简单，却蕴含着深刻的视觉原理：

x = X / (X + Y + Z)
y = Y / (X + Y + Z)

这个归一化过程正是CIE色度图的核心——它将三维色彩空间投影到二维平面，同时保留了人眼感知的色彩关系。

2. 构建光谱轨迹与色度图边界

CIE色度图的边界线（光谱轨迹）代表了纯单色光的颜色。要绘制这条曲线，我们需要计算各波长单色光的色度坐标：

def wavelength_to_xy(wavelength):
    """将特定波长转换为xy色度坐标"""
    idx = np.where(lambda_nm == wavelength)[0][0]
    X = x_bar[idx]
    Y = y_bar[idx]
    Z = z_bar[idx]
    sum_XYZ = X + Y + Z
    return X/sum_XYZ, Y/sum_XYZ

# 计算光谱轨迹
spectral_locus = []
for wl in range(360, 831, 5):  # 每5nm采样一次
    x, y = wavelength_to_xy(wl)
    spectral_locus.append((x, y))
spectral_locus = np.array(spectral_locus)

色度图内部的"马鞍形"区域包含了所有可见颜色，其中：

• 光谱色 ：位于边界线上的颜色，对应纯单色光
• 非光谱色 ：内部区域，需要通过混合不同波长的光才能产生
• 紫色线 ：连接380nm和780nm的直线，代表不存在的"纯紫色"光谱色

3. 可视化实现与Matplotlib技巧

使用Matplotlib绘制完整的CIE色度图需要一些技巧。以下是关键步骤的实现：

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch

fig, ax = plt.subplots(figsize=(10, 8))

# 绘制光谱轨迹
ax.plot(spectral_locus[:,0], spectral_locus[:,1], 'k-', lw=2)

# 绘制紫色线
purple_line = np.array([spectral_locus[0], spectral_locus[-1]])
ax.plot(purple_line[:,0], purple_line[:,1], 'k-', lw=2)

# 填充色度图区域
vertices = np.vstack([spectral_locus, purple_line[::-1]])
path = Path(vertices)
patch = PathPatch(path, facecolor='none', edgecolor='none')
ax.add_patch(patch)

# 创建网格和标签
ax.set_xlim(0, 0.8)
ax.set_ylim(0, 0.9)
ax.set_xlabel('x chromaticity')
ax.set_ylabel('y chromaticity')
ax.set_title('CIE 1931 Chromaticity Diagram')
ax.grid(True, alpha=0.3)

提示：要显示实际颜色，可以使用 imshow 叠加彩色渐变，或计算每个点的RGB近似值

4. 实际应用：色域可视化与色彩管理

理解色度图的最大价值在于解决实际开发中的色彩问题。我们可以用色度图来：

1. 比较不同设备的色域 ：

   def plot_gamut(ax, gamut_points, color, label):
       """在色度图上绘制设备色域"""
       path = Path(gamut_points)
       patch = PathPatch(path, facecolor=color, alpha=0.3, label=label)
       ax.add_patch(patch)
       ax.plot(gamut_points[:,0], gamut_points[:,1], color+'-', lw=1)
   
   # 示例：sRGB和Adobe RGB色域比较
   sRGB_points = np.array([[0.64,0.33], [0.30,0.60], [0.15,0.06]])
   adobeRGB_points = np.array([[0.64,0.33], [0.21,0.71], [0.15,0.06]])
   
   plot_gamut(ax, sRGB_points, 'r', 'sRGB')
   plot_gamut(ax, adobeRGB_points, 'b', 'Adobe RGB')
   ax.legend()
   

2. 颜色转换验证 ：

   def rgb_to_xy(rgb, color_space='sRGB'):
       """将RGB转换为xy色度坐标（简化版）"""
       if color_space == 'sRGB':
           # sRGB到XYZ的转换矩阵
           M = np.array([[0.4124, 0.3576, 0.1805],
                         [0.2126, 0.7152, 0.0722],
                         [0.0193, 0.1192, 0.9505]])
       # 其他色彩空间的矩阵...
       
       XYZ = np.dot(M, rgb)
       sum_XYZ = np.sum(XYZ)
       return XYZ[0]/sum_XYZ, XYZ[1]/sum_XYZ
   

3. 白点分析 ：

   # 常见标准白点
   white_points = {
       'D65': (0.3127, 0.3290),  # 标准日光
       'D50': (0.3457, 0.3585),  # 印刷行业常用
       'A': (0.4476, 0.4074)     # 白炽灯
   }
   
   for name, xy in white_points.items():
       ax.plot(xy[0], xy[1], 'ko')
       ax.annotate(name, xy=xy, xytext=(5,5), textcoords='offset points')
   

5. 高级话题：色差计算与色彩一致性

在实际开发中，我们经常需要量化两种颜色的差异。CIE色度图为此提供了基础：

def deltaE_CIE1976(xy1, xy2, Y1=1.0, Y2=1.0):
    """计算CIE1976 L*a*b*色差(简化版)"""
    # 首先转换到XYZ
    X1 = xy1[0] * Y1 / xy1[1]
    Z1 = (1 - xy1[0] - xy1[1]) * Y1 / xy1[1]
    
    X2 = xy2[0] * Y2 / xy2[1]
    Z2 = (1 - xy2[0] - xy2[1]) * Y2 / xy2[1]
    
    # 转换到Lab
    def xyz_to_lab(X, Y, Z, whitepoint):
        # 实现省略...
        return L, a, b
    
    L1, a1, b1 = xyz_to_lab(X1, Y1, Z1, white_points['D65'])
    L2, a2, b2 = xyz_to_lab(X2, Y2, Z2, white_points['D65'])
    
    return np.sqrt((L2-L1)**2 + (a2-a1)**2 + (b2-b1)**2)

注意：现代色彩科学使用更复杂的色差公式如CIEDE2000，但基本原理都源于CIE色度系统

6. 交互式可视化进阶

对于更深入的分析，我们可以创建交互式色度图：

from matplotlib.widgets import Slider

fig, ax = plt.subplots(figsize=(12, 9))
plt.subplots_adjust(bottom=0.25)

# 基本色度图绘制代码...

# 添加波长选择滑块
ax_wl = plt.axes([0.2, 0.1, 0.6, 0.03])
wl_slider = Slider(ax_wl, 'Wavelength (nm)', 360, 830, valinit=580, valstep=5)

def update(val):
    wl = wl_slider.val
    x, y = wavelength_to_xy(wl)
    point.set_data([x], [y])
    fig.canvas.draw_idle()

wl_slider.on_changed(update)

# 初始点
x, y = wavelength_to_xy(580)
point, = ax.plot(x, y, 'ro', markersize=10)

这种交互式可视化特别适合用于：

• 理解不同波长对应的颜色位置
• 演示颜色混合原理
• 分析特定颜色在色度图中的坐标

7. 性能优化与生产环境应用

当需要处理大量颜色数据时，纯Python实现可能效率不足。我们可以采用以下优化策略：

1. 向量化计算 ：

   def rgb_array_to_xy(rgb_array, color_space='sRGB'):
       """批量转换RGB到xy坐标"""
       if color_space == 'sRGB':
           M = np.array([[0.4124, 0.3576, 0.1805],
                         [0.2126, 0.7152, 0.0722],
                         [0.0193, 0.1192, 0.9505]])
       
       XYZ = np.dot(rgb_array, M.T)
       sum_XYZ = np.sum(XYZ, axis=1, keepdims=True)
       xy = XYZ[:, :2] / sum_XYZ
       return xy
   

2. 使用专业色彩科学库 ：

   import colour
   
   # 使用colour-science.org的专业实现
   xy = colour.models.RGB_to_xy(rgb_array, colour.models.RGB_COLOURSPACE_sRGB)
   

3. GPU加速 ：

   import cupy as cp
   
   def rgb_to_xy_gpu(rgb_array):
       """使用CuPy进行GPU加速计算"""
       rgb_gpu = cp.asarray(rgb_array)
       M = cp.array([[0.4124, 0.3576, 0.1805],
                     [0.2126, 0.7152, 0.0722],
                     [0.0193, 0.1192, 0.9505]])
       
       XYZ = cp.dot(rgb_gpu, M.T)
       sum_XYZ = cp.sum(XYZ, axis=1, keepdims=True)
       xy = XYZ[:, :2] / sum_XYZ
       return cp.asnumpy(xy)
   

在实际项目中，我发现将色度计算封装为单独的服务模块特别有用，这样可以在多个应用场景中复用。例如，在Web应用中可以通过REST API提供色彩转换服务，而在数据分析场景中则可以直接调用本地Python模块。