用MATLAB绘制CIE 1931色度图:从XYZ转换到色域叠加详解
发布时间:2026/9/13 12:59:47 作者:尧图编辑部 阅读量:1,286

简介这是一份面向色彩科学学习者、图像处理开发者以及MATLAB初学者的CIE色度图绘制脚本。资源包cie.rar内仅包含1个m文件压缩后大小约1KB代码量精简适合作为理解CIE标准配色与色度坐标的入门示例。目前已有704人学习下载。该脚本基于CIE 1931 XYZ色彩空间通过将X、Y、Z三刺激值转换为xy色度坐标在二维平面上绘制标准色度图并借助MATLAB绘图功能直观呈现色度与饱和度的关系。读者可直接运行或修改参数观察不同颜色在色度图上的位置变化也可在此基础上扩展交互查询、标注白点与光谱轨迹等功能。对于颜色匹配、显示校准、图像色彩分析等实际任务这份代码能帮助快速搭建可视化工具节省从头编写坐标转换与绘图流程的时间同时也有助于深入理解色度图背后的色彩科学原理。1. 拿到cie.m之后我第一件事是重画光谱轨迹下载过CIE色度图MATLAB代码的人大概都有同感压缩包很小打开就是一个几十行的cie.m脚本但直接运行往往和论文里的图对不上——边界发虚、颜色偏灰、白点位置模棱两可。我自己第一次跑这段代码也踩了这个坑后来把CIE 1931标准观察者数据、归一化公式、绘图管线逐行拆开核对才算真正把这张CIE色度图吃透。这篇文章不打算复述教科书结论而是把cie.m里涵盖的XYZ到xy坐标转换、光谱轨迹绘制、色域叠加三个环节拆开讲每段都附上能直接运行的MATLAB代码和参数调整建议。适合做图像处理、显示面板评估、LED照明或色彩管理的开发人员也适合研究生阶段需要自己绘制色度图做分析的读者。2. CIE 1931 XYZ到xy坐标转换色度坐标的计算逻辑2.1 从三刺激值到xy归一化不是可选项CIE 1931色彩空间的核心是三个线性光度量X、Y、Z分别对应人眼视网膜上三类视锥细胞在可见光谱范围内的加权响应。在cie.m这类MATLAB脚本里你不会直接对着光谱写积分而是预先加载一张380nm到780nm的色匹配函数表表里每一行是波长、x_bar、y_bar、z_bar四个值。拿到X、Y、Z之后第一步就是把三刺激值归一化成xy色度坐标这是色度图绘制的起点。公式很短但值得单独写出来x X / (X Y Z) y Y / (X Y Z)这里有一个隐藏风险如果光源数据来自非标准观察者文件或者积分区间太窄导致XYZ趋近于0归一化的除法会生成NaN绘图时对应点会直接消失马蹄形边界上就会留下缺口。我一般在归一化之前做一次分母保护常见做法是denom X Y Z; denom(denom 0) eps; % 用eps替换非正值防止NaN x X ./ denom; y Y ./ denom;这段代码的逻辑不复杂先求出三刺激值之和再用逻辑索引把非正数替换成eps最后做逐元素点除。参数上eps在双精度下约等于2.22e-16相比色度坐标的取值精度足够小既不会改变计算结果也不会触发除零警告。如果你处理的是单精度数组可以用realmin(single)替代eps原理一样。2.2 光谱轨迹与马蹄形边界是如何生成的光谱轨迹不是凭空“画”出来的而是逐波长“算”出来的。对每个波长λ把色匹配函数中对应行的xyz_bar归一化就得到该单色光的色度坐标(x(λ), y(λ))。把380nm到780nm所有波长按顺序连起来再用闭合曲线处理就形成经典的马蹄形边界。cie.m中常见的绘图顺序是先plot轮廓再fill填充也有人反过来先fill后plot区别仅在于边界线是否会被填充色盖住。波长间隔直接影响轨迹的平滑程度。标准CIE 1931 2度观察者数据通常是每5nm一个采样点380到780nm共81个点。如果你拿到的是10nm间隔的旧版本数据建议先做一次线性插值再绘图否则折线段会非常明显。插值不会改变色匹配函数的物理含义只是增加显示密度。用interp1即可完成注意把插值点范围保持在与原表一致。为了让验证有据可依我把几个关键波长的色度坐标列出来你可以直接对比自己计算的结果波长(nm)xy4500.15660.01775500.30160.69236000.62700.37296500.72600.2737看到这几个点大致吻合就说明XYZ数据源没有偏后续画出的马蹄形基本不会走样。2.3 白点坐标与参考基准的选择色度图只有轨迹边界还不够工程上需要白点作为中性参考。白点决定了色度学中“纯度”的距离基准也直接影响色域三角形在色度图中的相对位置。cie.m绘制完成后我习惯把D65白点单独标注出来。常用的标准白点坐标如下白点相关色温(K)xyD6565040.31270.3290D5050030.34570.3585A28560.44760.4074C67740.31010.3162这些数值来自CIE标准定义不是经验近似。D65是显示器、图像处理默认的日光白D50常见于印刷和色彩评审环境A对应白炽灯光源C是荧光灯时代的参考。在cie.m脚本里叠加任何色域之前先确认目标白点是哪一个再去选择对应的XYZ到RGB变换矩阵否则后续颜色映射会整体偏移。3. 用MATLAB绘制CIE色度图cie.m代码拆解与参数调整3.1 从标准观察者数据到光谱轨迹这里给出一份简化但仍可运行的cie.m实现。采用最常见的实现思路readtable读取CIE 1931标准观察者CSV归一化后用fill和plot描出马蹄形边界。你可以直接复制运行前提是当前目录下有格式正确的ciexyz31.csv。% 加载CIE 1931 2度标准观察者数据 % CSV格式第一列波长第二到第四列 x_bar y_bar z_bar data readtable(ciexyz31.csv); wavelength data.Wavelength; xBar data.x_bar; yBar data.y_bar; zBar data.z_bar; % 归一化到xy色度坐标注意是逐元素点除 denom xBar yBar zBar; x xBar ./ denom; y yBar ./ denom; % 先填充内部再描边界确保边界线清晰 figure(Color, w, Position, [100 100 720 640]); fill(x, y, [0.92 0.92 0.92], EdgeColor, none); hold on; plot(x, y, k-, LineWidth, 1.5); axis equal; grid on; xlabel(x (色度坐标)); ylabel(y (色度坐标)); title(CIE 1931 色度图);变量含义和绘图顺序值得展开说。readtable返回的是table对象后续用data.Wavelength取列比csvread这类旧函数直观得多。归一化那两行用的是点除MATLAB对矩阵运算非常严格写成/就是解方程程序会报错或给出错误结果。fill会先铺满闭合区域hold on确保后续plot叠加在同一坐标轴上最后的plot负责把黑色边界线重新描一遍否则fill的颜色会盖住边界。提示fill会自动把首尾两点连接成直线这条连线就是色度图中的紫线。plot只画实际光谱轨迹所以380nm和780nm两端不会自动连接两者之间的区域由fill负责闭合。3.2 给色度图铺一层彩色底图灰色填充图对阅读够用但要做报告或者拿给客户看最好在色度图内部铺上彩色。常见做法是把xy平面切分成网格每个点通过一个XYZ到RGB的近似矩阵映射成显示颜色再用imagesc输出。下面这段代码按sRGB D65标准来生成彩色底图% 生成300x300的xy网格 N 300; xv linspace(0, 0.8, N); yv linspace(0, 0.9, N); [Xg, Yg] meshgrid(xv, yv); Zg 1 - Xg - Yg; % 只有x、y、z都非负且和为1的点才落在马蹄形内 valid (Xg 0) (Yg 0) (Zg 0) (Zg 1); % sRGB D65的XYZ到线性RGB变换矩阵 M [ 3.2406 -1.5372 -0.4986; -0.9689 1.8758 0.0415; 0.0557 -0.2040 1.0570]; RGB zeros(N, N, 3); for c 1:3 RGB(:, :, c) M(c,1) * Xg M(c,2) * Yg M(c,3) * Zg; end RGB max(0, min(1, RGB)); % 裁剪超出[0,1]的显示值 % 非马蹄形区域统一填白色逐通道处理 img RGB; mask ~valid; for c 1:3 channel img(:, :, c); channel(mask) 0.95; img(:, :, c) channel; end % 可视化并叠加光谱轨迹 imagesc(xv, yv, img); axis xy; hold on; plot(x, y, k-, LineWidth, 2);这段代码有两个关键点。第一valid矩阵的构造条件是x、y、z都必须非负且xyz等于1其中z由1-x-y推得z非负意味着网格点必须落在线段xy≤1以内。超出有效区域的点不能显示颜色否则色度图会出现一大片没有物理意义的色块。第二M矩阵是sRGB D65下的标准线性变换它假设输入XYZ以D65白点归一化。如果你的数据基于D50画出来的颜色会整体偏暖必须先做Bradford色适应转换之后再去套这个矩阵。参数N300是底图分辨率数值越大颜色过渡越细腻但计算量按平方增长。普通笔记本跑300没有问题调试阶段建议先用150确认效果后再调高。M矩阵的每个系数都有明确物理含义第一行是XYZ到R的权重第二行到G第三行到B。输出的是线性RGB严格来说还需要再做一次sRGB gamma编码才能得到最终显示值但这里省略了gamma步骤因为底图仅用于示意肉眼对颜色细节的敏感度有限。3.3 坐标范围约束、波长标注与调试技巧色度图的坐标范围如果不固定每次重绘马蹄形的比例都会变叠图对比时会产生误导。常见做法是执行xlim([0 0.8])和ylim([0 0.9])两侧留白足够容纳标注文字。关键波长可以用text直接标记% 在光谱轨迹上标注关键波长 keyWavelengths [450 500 530 550 600 650 700]; for k 1:length(keyWavelengths) idx find(abs(wavelength - keyWavelengths(k)) 2.5); if ~isempty(idx) text(x(idx(1)) 0.01, y(idx(1)) 0.01, ... [num2str(keyWavelengths(k)) nm], FontSize, 8); end endfind配合abs窗口2.5nm是为了兼容波长表间隔不是严格5nm的情况。text的前两个参数是位置偏移量0.01让文字离开曲线边界。如果你希望图面更干净可以把偏移量加大到0.02但注意不要超出坐标范围。调试阶段常用的一组参数我整理在下面参数推荐值说明N300底图网格分辨率调试时用150LineWidth1.5 ~ 2光谱轨迹线宽导出图片前调大xlim / ylim[0 0.8] / [0 0.9]固定坐标范围便于叠图FontSize8 ~ 10波长标注字号4. 在cie.m基础上叠加色域三角形、白点与普朗克轨迹4.1 sRGB色域三角形与覆盖率计算色度图画完之后最常做的一件事就是把设备的色域三角形叠上去。sRGB色域在CIE 1931色度图上的三个顶点是固定的红色(0.6400, 0.3300)、绿色(0.3000, 0.6000)、蓝色(0.1500, 0.0600)。把三个点按顺序连成一个闭合三角形就能直观看到设备色彩范围和可见光谱的关系。cie.m脚本里一般没有现成的色域叠加函数需要自己用plot和fill补上。% 定义sRGB色域三个顶点最后一个点回到起点闭合图形 sRGB_Gamut [0.6400 0.3300; 0.3000 0.6000; 0.1500 0.0600; 0.6400 0.3300]; plot(sRGB_Gamut(:, 1), sRGB_Gamut(:, 2), r-, LineWidth, 1.8); hold on; fill(sRGB_Gamut(:, 1), sRGB_Gamut(:, 2), [1 0.9 0.9], ... EdgeColor, none, FaceAlpha, 0.3); % 计算sRGB在可见光谱范围内的近似覆盖率 inGamut inpolygon(Xg, Yg, sRGB_Gamut(:, 1), sRGB_Gamut(:, 2)); coverage sum(inGamut(:)) / sum(valid(:)); fprintf(sRGB在CIE1931色度图中的近似覆盖率: %.2f%%\n, coverage * 100);inpolygon是MATLAB自带的点集包含判定函数传入多边形顶点坐标后自动判断网格点是否落在内部。覆盖率计算中分子是sRGB三角形内的网格点数分母是马蹄形有效区域内的网格点数比值就是近似面积覆盖率。这个值依赖底图分辨率N但N超过200后变化已经小于0.1%用来横向比较不同色域标准完全够用。需要注意valid矩阵必须来自3.2节如果你跳过了底图生成就需要单独用inpolygon把马蹄形内部区域先算出来。4.2 不同色域标准的一次性对比查看覆盖率只是第一步实际工作中更多是在一张图上对比多种色域。sRGB、Adobe RGB、DCI-P3是显示行业最常见的三种基准它们的顶点坐标差别很大特别是绿色顶点色域标准红(x, y)绿(x, y)蓝(x, y)白点sRGB(0.640, 0.330)(0.300, 0.600)(0.150, 0.060)D65Adobe RGB(0.640, 0.330)(0.210, 0.710)(0.150, 0.060)D65DCI-P3(0.680, 0.320)(0.265, 0.690)(0.150, 0.060)D65Adobe RGB的绿色顶点y坐标达到0.710明显比sRGB的0.600更靠近光谱轨迹上缘这就是为什么它在绿色和青色区域能覆盖更多颜色。DCI-P3则把红色顶点外移到x0.680在红色区域表现更强。实际对比时可以写一个循环把多个色域依次画出来用不同颜色和线型区分。注意每画一个三角形之前都要执行hold on否则前一个色域会被清掉。4.3 叠加D65白点、普朗克轨迹与显示器色温判断白点和色温线叠加是色度图进阶用法里频率最高的需求主要用于判断显示器的色温预设是否偏色。在cie.m脚本末尾追加下面这段代码% 标注D65白点 d65 [0.3127 0.3290]; plot(d65(1), d65(2), ks, MarkerSize, 8, MarkerFaceColor, w); text(d65(1) 0.01, d65(2) 0.01, D65, FontSize, 10); % 绘制普朗克轨迹约1000K到10000K T linspace(1000, 10000, 91); xyPlanck zeros(length(T), 2); for i 1:length(T) xyPlanck(i, :) planckianLocus(T(i)); % 需要自己实现的函数 end plot(xyPlanck(:, 1), xyPlanck(:, 2), b-, LineWidth, 1.2);planckianLocus是cie.m里通常需要补全的函数MATLAB没有内置。实现思路是根据普朗克黑体辐射公式计算色温T对应的光谱功率分布再与CIE色匹配函数做内积得到XYZ最后归一化为xy坐标。如果不想自己写积分可以直接用CIE公开的黑体轨迹色度表通过interp1插值到任意色温误差在可接受范围内。实际项目里把一台显示器的实测白点坐标放到色度图上看它偏离普朗克轨迹的距离和方向。偏离轨迹上方偏品红下方偏绿左右偏移则反映色温误差。一般D65标准的显示器白点偏差在xy各0.003以内人眼基本无感超过0.01就很明显了。这个判断方法直接来自CIE色度学不会因为显示器厂商的软件校准而改变。5. 验证CIE色度图绘制是否正确的五个检查点5.1 用已知坐标校验边界和关键点调试色度图画得好不好最直接的办法是对数值而不是对感觉。取三个已知的色度坐标比如sRGB的红色顶点(0.6400, 0.3300)在图中查询对应位置的颜色是否符合预期。更可靠的是检查550nm边界点坐标理论值约为(0.3016, 0.6923)这是光谱轨迹最接近顶部的点如果你算出来的数据偏差超过0.005多半是色匹配函数表加载错误。5.2 检查光谱轨迹两端点相对位置光谱轨迹两个端点都在横轴附近380nm端点落在左下角紫色区域780nm端点收拢到红光末端一侧。如果画完图发现两端位置明显偏离优先排查CSV末尾是否有空行或被读成NaN的列常见原因是readtable把某些数值列当成了文本导致x_bar、y_bar无法参与归一化运算。5.3 用gamma修正让底图颜色更接近人眼感知彩色底图生成时省略了gamma编码步骤直接在线性RGB空间显示画面通常会偏暗。想要更接近显示器实际观感可以对RGB做一次标准sRGB gamma变换小于等于0.0031308的按12.92倍线性放大否则按1.055倍开1/2.4次方再减去0.055。这一步可以在imagesc之前完成代码量不超过五行但底图会明显更通透。5.4 把图导出成可复用的高质量结果报告或者论文插图建议用exportgraphics导出分辨率参数设300dpi格式选PNG或TIFF。如果直接在figure窗口截图像素密度往往不够印刷出来会虚。需要矢量化输出时可以用print -depsc生成EPS再交给LaTeX编译线宽不会因为缩放而改变。本文还有配套的精品资源点击获取