用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. 实际应用:色域可视化与色彩管理
理解色度图的最大价值在于解决实际开发中的色彩问题。我们可以用色度图来:
-
比较不同设备的色域 :
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() -
颜色转换验证 :
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 -
白点分析 :
# 常见标准白点 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实现可能效率不足。我们可以采用以下优化策略:
-
向量化计算 :
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 -
使用专业色彩科学库 :
import colour # 使用colour-science.org的专业实现 xy = colour.models.RGB_to_xy(rgb_array, colour.models.RGB_COLOURSPACE_sRGB) -
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模块。
1394

被折叠的 条评论
为什么被折叠?



