视频编码自学教程 02：颜色与色彩空间——从 RGB 到 Y'CbCr 4:2:0用手机拍一张照片，屏幕上看到的是红、绿

系列导航：目录｜上一章：视频是什么｜下一章：视频压缩全景（编写中）

本章目标：读完本章，你将能够

说清楚为什么三个数就能表示一种颜色，读懂 CIE 色度图，比较 BT.709、DCI-P3、BT.2020 色域的大小；

由三原色和白点推导出 RGB → XYZ 矩阵，知道亮度公式里 0.2126、0.7152、0.0722 是怎么来的；

理解 gamma（传递函数）为什么必不可少，用 ΔL* 定量比较线性编码和 gamma 编码，了解 HDR 的 PQ 和 HLG；

熟练地在 R'G'B' 和 Y'CbCr 之间换算，知道视频里的黑为什么是 16、白为什么是 235；

理解 4:2:2、4:2:0 的含义和色度样本的位置，算出色度二次采样损失了多少；

认出范围用错、矩阵用错这两类最常见的颜色事故，读写 I420、NV12 等原始 YUV 文件。

前置知识：第 01 章的采样、量化和 PSNR；矩阵乘法与 3×3 线性方程组（会用 numpy 解即可）；对数和幂函数。

用手机拍一张照片，屏幕上看到的是红、绿、蓝三种小发光点组成的图像；把它存成 H.264 视频，存下来的却不是 R、G、B，而是一个“亮度”加两个“色差”，而且色差只有一半甚至四分之一的分辨率。播放时再换算回 RGB 显示出来。

为什么要绕这么大一个圈？简单地说，是为了迁就人眼：人眼对亮度的细节很敏感，对颜色的细节却很迟钝；人眼对亮度的感受是对数式的，而不是线性的。视频编码在真正开始“压缩”之前，就已经用这两点省下了一半的数据（4:4:4 的 24 bit/像素降到 4:2:0 的 12 bit/像素）。

这一章把颜色从物理讲到文件：光和眼睛、色度图和色域、gamma、Y'CbCr、色度二次采样，最后落到硬盘上的 YUV 文件。几乎每一个视频工程师都踩过颜色的坑：画面发灰、偏色、黑色发紫……读完本章，你会知道这些问题从哪里来。

一、颜色是什么：光、眼睛与三个数

1. 光谱与三种视锥细胞

可见光是波长约 380～780 nm 的电磁波。一束光的物理性质由它的光谱 $S(\lambda)$ 完全描述：每个波长上有多少能量。理论上这是一个函数，有无穷多个自由度。

但人眼并不测量完整的光谱。视网膜上负责白天视觉的视锥细胞只有三种，分别对长、中、短波长最敏感（L、M、S 视锥）。每种视锥细胞只输出一个数：它对整个光谱的加权响应

L = \int S(\lambda)\, l(\lambda)\,\mathrm{d}\lambda,\qquad M = \int S(\lambda)\, m(\lambda)\,\mathrm{d}\lambda,\qquad S_{\text{视锥}} = \int S(\lambda)\, s(\lambda)\,\mathrm{d}\lambda

无穷维的光谱经过眼睛，被“投影”成了三个数。这就是三色原理（trichromacy）：人看到的颜色只由这三个数决定。

2. 同色异谱：为什么三个数就够

既然只剩三个数，那么两束光谱完全不同的光，只要这三个积分相等，看起来就是同一种颜色。这叫同色异谱（metamerism）。

同色异谱是所有彩色显示技术的基础。屏幕不可能复现自然界里任意的光谱，它只有红、绿、蓝三种发光点。但它不需要复现光谱，只要调节三种发光点的强度，让眼睛的三个响应和原来的光一样就行了。所以颜色本质上是一个三维量，RGB、XYZ、Y'CbCr 只是这个三维空间里的不同坐标系。

3. CIE XYZ 与配色函数

1931 年，国际照明委员会（CIE）根据大量配色实验，定义了一套标准的“虚拟眼睛”：CIE 1931 标准观察者。它用三条配色函数 $\bar x(\lambda)$ 、 $\bar y(\lambda)$ 、 $\bar z(\lambda)$ 代替真实的视锥响应：

X = \int S(\lambda)\,\bar x(\lambda)\,\mathrm{d}\lambda,\qquad Y = \int S(\lambda)\,\bar y(\lambda)\,\mathrm{d}\lambda,\qquad Z = \int S(\lambda)\,\bar z(\lambda)\,\mathrm{d}\lambda

图 1：CIE 1931 标准观察者的三条配色函数。CIE 公布的是数值表格，这里用 Wyman 等人 2013 年提出的多段高斯近似公式计算（附录代码第 1 节）。下方色带是各波长单色光的大致颜色。

XYZ 不是真实的视锥响应，而是它们的一个线性组合。这个组合是特意挑选的，有两个好性质：

三条曲线处处非负，所以任何真实光的 X、Y、Z 都不会是负数；
$\bar y(\lambda)$ 恰好等于人眼的明视觉光谱光效率函数，所以 $Y$ 就是物理意义上的亮度（luminance）。 $\bar y$ 在约 555 nm 的黄绿光处最高（近似公式给出 554 nm），说明人眼对黄绿光最敏感，对蓝光和深红光很不敏感。

记住 $Y$ 这个字母，它会贯穿整章，最后变成 Y'CbCr 里的那个 Y'。

二、色度图与色域

1. 色度坐标

XYZ 三个数里混着两类信息：多亮，以及什么颜色。把亮度除掉，就只剩“什么颜色”，叫色度（chromaticity）：

x = \frac{X}{X+Y+Z},\qquad y = \frac{Y}{X+Y+Z},\qquad z = 1 - x - y

只需要 $(x,\ y)$ 两个数。反过来，已知色度 $(x,\ y)$ 和亮度 $Y$ ，就能还原

X = \frac{x}{y}\,Y,\qquad Z = \frac{1-x-y}{y}\,Y

把每个波长的单色光算出 $(x,\ y)$ 画在平面上，就得到一条马蹄形曲线，叫光谱轨迹。比如近似公式算出 520 nm 的绿光在 $(0.0813,\ 0.8197)$ ，620 nm 的红光在 $(0.6955,\ 0.3045)$ 。连接两端的直线叫紫线，上面是红和蓝的混合色，光谱里没有单一波长对应它们。所有真实颜色都在马蹄形之内。

图 2：CIE 1931 xy 色度图和三种常见色域。光谱轨迹用近似公式计算，只画了 420～650 nm；更短或更长的波长色度坐标几乎不再变化，挤在马蹄形的两个端点上。图中的颜色只是示意，屏幕显示不出 BT.709 三角形以外的颜色。

2. 三原色与白点

一台显示器的三种发光点各有一个色度坐标，叫三原色（primaries）。三种原色按不同强度混合，能得到的所有颜色就是以三个原色为顶点的三角形，叫这台设备的色域（gamut）。三角形以外的颜色，无论怎么调都显示不出来。

R = G = B = 1（三种原色全开）时混合出的颜色叫白点。视频和电脑显示几乎都用 D65，色度 $(0.3127,\ 0.3290)$ ，大致相当于正午的日光。

标准	红 $(x,\ y)$	绿 $(x,\ y)$	蓝 $(x,\ y)$	用途
BT.709 / sRGB	(0.640, 0.330)	(0.300, 0.600)	(0.150, 0.060)	高清电视、网络视频、绝大多数电脑屏幕
DCI-P3	(0.680, 0.320)	(0.265, 0.690)	(0.150, 0.060)	数字电影；苹果等设备的 Display P3 用同样的原色配 D65 白点
BT.2020	(0.708, 0.292)	(0.170, 0.797)	(0.131, 0.046)	超高清电视、HDR

（数字电影院的 DCI-P3 原本用另一个偏绿的白点，这里和附录代码统一用 D65。）

3. 色域有多大

用三角形面积比较色域（附录代码第 3 节）：

色域	xy 面积	相当于 BT.2020 的	u'v' 面积	相当于 BT.2020 的
BT.709	0.1120	52.9%	0.0649	58.0%
DCI-P3	0.1520	71.7%	0.0815	72.9%
BT.2020	0.2119	100%	0.1118	100%

xy 色度图有一个缺点：在绿色区域，坐标差很远的两种颜色看起来可能差不多；在蓝紫区域，坐标差一点点颜色就变了很多。CIE 1976 年给出了一个改进的 $u'v'$ 坐标，更接近“距离相等，看起来差别也相等”：

u' = \frac{4x}{-2x + 12y + 3},\qquad v' = \frac{9y}{-2x + 12y + 3}

按 $u'v'$ 面积算，BT.709 约占 BT.2020 的 58%。高清电视用了二十多年的 BT.709，只能显示 BT.2020 色域的一半多一点。BT.2020 的三个原色都在光谱轨迹上，是纯的单色光，所以它的三角形几乎贴着马蹄形的边。

三、RGB 与 XYZ 的换算：亮度系数从哪里来

1. 推导

线性 RGB（R、G、B 与发光强度成正比）和 XYZ 之间是一个 3×3 线性变换：

\begin{bmatrix} X \\ Y \\ Z \end{bmatrix} = \mathbf{M} \begin{bmatrix} R \\ G \\ B \end{bmatrix}

$\mathbf{M}$ 完全由三原色的色度和白点决定，推导分两步。

第一步：方向。 只开红色（ $R = 1$ ， $G = B = 0$ ）时，输出的 XYZ 就是 $\mathbf{M}$ 的第一列。它的色度必须是红原色 $(x_r,\ y_r)$ ，所以第一列一定正比于亮度为 1 时的 XYZ：

\mathbf{p}_r = \left(\frac{x_r}{y_r},\ 1,\ \frac{1 - x_r - y_r}{y_r}\right)^\top

绿、蓝两列同理。于是 $\mathbf{M} = \begin{bmatrix} S_r\mathbf{p}_r & S_g\mathbf{p}_g & S_b\mathbf{p}_b \end{bmatrix}$ ，只剩三个缩放系数 $S_r$ 、 $S_g$ 、 $S_b$ 未知。

第二步：用白点定缩放。 三原色全开（ $R = G = B = 1$ ）要得到亮度为 1 的白点 $\mathbf{w}$ ：

S_r\,\mathbf{p}_r + S_g\,\mathbf{p}_g + S_b\,\mathbf{p}_b = \mathbf{w} \quad\Longrightarrow\quad \begin{bmatrix} \mathbf{p}_r & \mathbf{p}_g & \mathbf{p}_b \end{bmatrix} \begin{bmatrix} S_r \\ S_g \\ S_b \end{bmatrix} = \mathbf{w}

解这个 3×3 线性方程组就得到 $\mathbf{M}$ 。代入 BT.709 的原色和 D65：

\mathbf{M}_{709} = \begin{bmatrix} 0.4124 & 0.3576 & 0.1805\\ 0.2126 & 0.7152 & 0.0722\\ 0.0193 & 0.1192 & 0.9505 \end{bmatrix}

2. 三套亮度系数

看 $\mathbf{M}$ 的第二行： $Y = 0.2126\,R + 0.7152\,G + 0.0722\,B$ 。因为每个 $\mathbf{p}$ 的中间分量都是 1，第二行恰好就是 $(S_r,\ S_g,\ S_b)$ ，也就是三原色全开组成白光时，各自贡献的亮度。

这就是视频里亮度公式的来历。它们不是凭经验拍脑袋定的，而是由原色和白点算出来的。换一套原色，系数就跟着变（附录代码第 2 节）：

原色 + 白点	$K_R$	$K_G$	$K_B$	用在
NTSC 1953 + C 光源	0.2989	0.5866	0.1144	BT.601（标清）、JPEG
BT.709 + D65	0.2126	0.7152	0.0722	BT.709（高清）
BT.2020 + D65	0.2627	0.6780	0.0593	BT.2020（超高清）

第 01 章用的 $0.299,\ 0.587,\ 0.114$ 来自 1953 年 NTSC 彩色电视的原色。有意思的是，后来的标清电视实际上早已换用了别的原色，但 BT.601 和 JPEG 一直沿用这组老系数。这是历史包袱，也是第八节“矩阵用错”问题的根源。

3. 色域之间的转换

有了两个色彩空间各自的 $\mathbf{M}$ ，从 BT.709 转到 BT.2020 就是 $\mathbf{M}_{2020}^{-1}\mathbf{M}_{709}$ ：

\begin{bmatrix} R \\ G \\ B \end{bmatrix}_{2020} = \begin{bmatrix} 0.6274 & 0.3293 & 0.0433\\ 0.0691 & 0.9195 & 0.0114\\ 0.0164 & 0.0880 & 0.8956 \end{bmatrix} \begin{bmatrix} R \\ G \\ B \end{bmatrix}_{709}

这和 ITU-R BT.2087 标准里给出的矩阵四位小数完全一致，可以作为推导正确的验证。几个要点：

这个变换作用在线性光上，所以 gamma 编码的信号必须先解码成线性值，转换完再重新编码；
BT.709 的纯红 $(1, 0, 0)$ 在 BT.2020 里是 $(0.6274,\ 0.0691,\ 0.0164)$ 。小色域的颜色总能在大色域里找到；
反过来，BT.2020 的纯绿在 BT.709 里是 $(-0.5876,\ 1.1329,\ -0.1006)$ ，出现了负数和大于 1 的数。显示器发不出“负的红光”，这个颜色超出了 BT.709 色域。把它塞进小色域叫色域映射（gamut mapping），最简单的办法是截断到 $[0,\ 1]$ ，但会造成颜色失真，更好的方法是第 27 章讲 HDR 蓝光时会遇到的色调映射。

四、传递函数：为什么要做 gamma

1. 人眼对亮度的感受是非线性的

把一盏灯的亮度从 1 加倍到 2，和从 100 加倍到 200，人感觉到的“变亮了多少”差不多。也就是说，人眼对亮度的相对变化敏感，对绝对变化不敏感。CIE 用明度（lightness） $L^*$ 来描述这种感受：

L^* = \begin{cases} 116\,Y^{1/3} - 16, & Y > (6/29)^3 \approx 0.008856\\[2pt] 903.3\,Y, & \text{否则} \end{cases}

其中 $Y$ 是相对于白色的亮度（0～1）， $L^*$ 从 0（黑）到 100（白）。 $L^*$ 大致是亮度的立方根：亮度只有白色 18% 的灰， $L^*$ 已经接近 50，看起来是“中灰”。两种颜色的 $L^*$ 相差约 1 时，人眼大约刚好能分辨出来。

2. OETF 与 EOTF

既然人眼的感受是非线性的，存储亮度时就不应该线性地分配码值，而应该在暗部分得更细、亮部分得更粗。做法是在摄像机里先把线性光 $L$ 经过一条压缩曲线变成编码值 $V$ ，叫光电转换函数（Opto-Electronic Transfer Function，OETF）。BT.709 的 OETF 是

V = \begin{cases} 4.5\,L, & L < 0.018\\ 1.099\,L^{0.45} - 0.099, & L \ge 0.018 \end{cases}

显示端用电光转换函数（EOTF）把 $V$ 还原成发光强度。BT.1886 规定高清显示器的 EOTF 近似为 $L = V^{2.4}$ 。电脑图像常用的 sRGB 曲线和它形状相近。

图 3：(a) 几条传递函数和明度 $L^*/100$ ，gamma 曲线和人眼的感受曲线形状很接近；(b) 8 位编码下相邻两级码值的明度差 ΔL。*

这条曲线俗称 gamma，因为它近似一个幂函数 $V = L^{1/\gamma}$ 。从历史上看，它是因为 CRT 显像管的电压和亮度之间天然就是约 2.5 次方的关系，摄像机必须预先做反向校正。但它被保留到今天，并不只是为了兼容老电视：碰巧，CRT 的非线性和人眼的感受差不多正好相反，gamma 编码恰好让码值在“人眼感受”上分布得比较均匀。

还要注意，摄像机的 OETF（约 0.5 次方）和显示器的 EOTF（2.4 次方）并不互逆，端到端的总效果约是 1.2 次方。这是有意为之的：在昏暗的客厅里看电视，适度提高对比度看起来更接近现场。

3. 8 位线性编码不够用

用 ΔL* 可以定量地看 gamma 的作用（附录代码第 4 节）。把 0～1 的亮度分成 256 级：

编码方式	第一级的 ΔL*	最大 ΔL*	ΔL* > 1 的台阶数	$L^* < 50$ 的暗部分到的码值
8 位线性	3.54	3.54	15	47 个
8 位 sRGB	0.27	0.51	0	119 个

线性编码把大部分码值浪费在人眼分辨不出差别的亮部，暗部却只有 47 个码值，最暗的十几个台阶每一步都超过了“刚好能察觉”的差别。sRGB 编码把一半的码值分给了暗部，每一步都在 ΔL* = 0.51 以内。附录代码还算出：线性编码要用到 11 位，最大台阶才能和 8 位 sRGB 一样小。gamma 编码相当于免费多出了 3 位。

图 4：(a) 从黑到中灰的渐变，8 位线性编码在暗部出现明显的台阶，8 位 sRGB 编码看不出台阶；(b) HDR 的 PQ 曲线。

所以，视频编码器拿到的数据几乎都是 gamma 编码过的，编码器只是把它们当作一堆整数来压缩，并不关心它们代表多亮。只有在需要做物理意义上的计算时，比如色域转换、混合两种颜色、高质量缩放，才需要先解码成线性光。

4. 撇号的含义：luma 与 luminance

从这里开始，字母上的撇号很重要：

$R$ 、 $G$ 、 $B$ 、 $Y$ 表示线性的量， $Y$ 是真正的亮度（luminance）；
$R'$ 、 $G'$ 、 $B'$ 表示经过 gamma 编码的值，用它们加权求和得到的

Y' = K_R R' + K_G G' + K_B B'

叫 luma，中文里也常叫“亮度”，但它并不等于真正亮度的 gamma 编码值。对灰色两者相等，对饱和的颜色会有偏差。这样做是工程上的妥协：直接用 $R'G'B'$ 算，电路简单，代价是一小部分亮度信息“漏”进了色度通道。BT.2020 定义了一种“恒定亮度”的变体来纠正它，但实际几乎没人用。本教程沿用行业习惯：说“亮度”时一般指 $Y'$ 。

5. HDR：PQ 与 HLG

传统的 BT.709 是为最高约 100 cd/m²（尼特）的显示器设计的。HDR（高动态范围）视频要表示从 0.001 到上万尼特的亮度，传递函数也必须换掉。BT.2100 标准规定了两种：

PQ（Perceptual Quantizer，SMPTE ST 2084）：编码的是绝对亮度，最高 10000 cd/m²。曲线按人眼在不同亮度下的分辨能力设计，10 位就足够。附录代码算出，100 cd/m² 编码成 0.508，1000 cd/m² 编码成 0.752，也就是说 SDR 的全部亮度范围只用了 PQ 一半左右的码值；
HLG（Hybrid Log-Gamma，BBC 和 NHK 提出）：编码的是相对亮度，暗部用平方根曲线（与传统 gamma 相近），亮部换成对数曲线，分界点处的编码值恰好是 0.5。它的好处是不认识 HDR 的老电视也能大致正常显示。

HDR10、杜比视界用 PQ，广播电视多用 HLG。后面讲 H.265 和 UHD 蓝光时还会遇到它们。

五、Y'CbCr：一个亮度加两个色差

1. 定义

有了 $Y'$ ，颜色信息就用两个色差来表示： $B' - Y'$ 和 $R' - Y'$ 。灰色的 $R' = G' = B'$ ，两个色差都是 0。为了让它们的取值范围恰好是 $[-0.5,\ 0.5]$ ，再各除以一个常数：

P_B = \frac{B' - Y'}{2\,(1 - K_B)},\qquad P_R = \frac{R' - Y'}{2\,(1 - K_R)}

为什么是这两个分母？以 $P_B$ 为例， $B' - Y' = (1 - K_B)B' - K_R R' - K_G G'$ 。 $R'$ 、 $G'$ 、 $B'$ 都在 $[0,\ 1]$ 内，所以 $B' - Y'$ 最大是 $1 - K_B$ （纯蓝），最小是 $-(1 - K_B)$ （黄色， $R' = G' = 1$ ， $B' = 0$ ），除以 $2(1 - K_B)$ 后正好落在 $[-0.5,\ 0.5]$ 内。

$P_B$ 、 $P_R$ 是模拟或浮点形式。量化成整数后叫 $C_B$ 、 $C_R$ （下一节），通常写作 Cb、Cr。整个色彩空间叫 Y'CbCr。

2. 矩阵与逆矩阵

代入 BT.709 的系数，写成矩阵：

\begin{bmatrix} Y' \\ P_B \\ P_R \end{bmatrix} = \begin{bmatrix} 0.2126 & 0.7152 & 0.0722\\ -0.1146 & -0.3854 & 0.5\\ 0.5 & -0.4542 & -0.0458 \end{bmatrix} \begin{bmatrix} R' \\ G' \\ B' \end{bmatrix}

每一行的和分别是 1、0、0：白色映射到 $(1,\ 0,\ 0)$ ，所有灰色的色差都是 0。逆变换形式很简单：

R' = Y' + 1.5748\,P_R,\qquad G' = Y' - 0.1873\,P_B - 0.4681\,P_R,\qquad B' = Y' + 1.8556\,P_B

其中 $1.5748 = 2(1 - K_R)$ ， $1.8556 = 2(1 - K_B)$ 。BT.601 和 BT.2020 的矩阵只是换了 $K_R$ 、 $K_B$ ，附录代码第 5 节把三套矩阵都打印出来了。

3. 能量集中到亮度

图 5：上排是 R'、G'、B' 三个通道，下排是 Y'、Cb、Cr。R'G'B' 三张图几乎一样；Y' 包含了全部细节，两个色差通道大片是中性的灰，只在有颜色的地方有内容（色差通道的显示对比度放大了 2 倍）。

图 5 一眼就能看出 Y'CbCr 的好处：R'G'B' 三个通道高度相似，每一张都带着完整的轮廓和细节，存三份是很大的浪费；换成 Y'CbCr 后，细节几乎都集中到了 Y'。附录代码第 6 节定量地算了出来：

	通道两两之间的相关系数	三个通道各占总方差的比例
R'G'B'	R'G' 0.84，R'B' 0.73，G'B' 0.96	36.1%，31.5%，32.5%
Y'CbCr	Y'Cb −0.08，Y'Cr −0.05，CbCr −0.77	88.6%，3.1%，8.3%

图 6：(a) G' 和 B' 几乎落在一条直线上；(b) Y' 和 Cr 之间看不出线性关系；(c) 三个通道的方差占比，Y'CbCr 和主成分分析都把大约 89% 的能量集中到了第一个通道。

这和第 01 章的“相邻像素很像”是同一个道理，只是换成了“同一个像素的三个通道很像”。消除这种相关性叫去相关，是一切变换编码的出发点。

两个细节值得说一下：

Cb 和 Cr 之间仍然相关（−0.77）。这张照片里大面积的橙色宇航服和肤色，Cb 都偏负、Cr 都偏正。Y'CbCr 是固定的矩阵，不会根据图像调整，所以不可能把每张图都去相关得很彻底。
主成分分析（PCA）是“针对这张图最优”的去相关变换：它找到的第一个方向几乎是 R'、G'、B' 的平均 $(0.33,\ 0.34,\ 0.33)$ ，占总方差的 89.3%，只比 Y' 的 88.6% 多一点点。可见固定的 Y'CbCr 已经很接近最优，而且它的 Y' 对应人眼的亮度感受，这一点 PCA 做不到。第 07 章讲 DCT 时，还会遇到 PCA 的另一个名字：KLT（Karhunen–Loève 变换）。

图 6(b) 还提醒我们：相关系数接近 0 只说明没有线性关系，Y' 和 Cr 之间显然还有别的结构。H.266 的“跨分量线性模型”（CCLM）就利用了亮度和色度之间的局部关系来预测色度（第 25 章）。

4. “YUV”这个名字

在软件和日常交流里，Y'CbCr 几乎总是被叫作 YUV，FFmpeg 的像素格式名也是 yuv420p。严格地说，YUV 是 PAL 模拟电视的色彩空间， $U = 0.492\,(B' - Y')$ ， $V = 0.877\,(R' - Y')$ ，缩放系数和 Y'CbCr 不同。但在数字视频的语境里，“YUV”基本就是 Y'CbCr 的同义词，本教程也会混用：说“YUV 文件”时，指的就是存着 Y'CbCr 整数的文件。

六、量化与范围：为什么黑是 16

$Y'$ 在 $[0,\ 1]$ ， $P_B$ 、 $P_R$ 在 $[-0.5,\ 0.5]$ ，存储前要量化成整数。视频里最常用的是有限范围（limited range，也叫 TV range、studio range）：

Y = \operatorname{round}\!\big(2^{N-8}\,(219\,Y' + 16)\big),\qquad C = \operatorname{round}\!\big(2^{N-8}\,(224\,P + 128)\big)

$N$ 是位深。8 位时，黑是 16，白是 235，色差从 16 到 240，128 表示没有颜色；10 位时乘以 4，黑是 64，白是 940。

为什么不用满 0～255？主要有两个原因：

留出余量：视频在制作和传输过程中要经过各种滤波，锐利的边缘处会产生过冲和下冲（第 01 章波带片的例子里见过）。有限范围在黑以下和白以上留了余量，过冲的值不会立刻被截断；
保留特殊码值：在演播室使用的串行数字接口（SDI，BT.656 规定）里，8 位的 0 和 255 被保留下来，作为每一行开始和结束的定时标记，不能出现在图像数据里。

另一种是全范围（full range，也叫 PC range）： $Y = 255\,Y'$ ， $C = 255\,P + 128$ 。JPEG 图片用全范围，电脑屏幕截图和游戏录屏也常用全范围。

100% 彩条按 BT.709、8 位有限范围编码（附录代码第 5 节）：

颜色	Y	Cb	Cr
白	235	128	128
黄	219	16	138
青	188	154	16
绿	173	42	26
品红	78	214	230
红	63	102	240
蓝	32	240	118
黑	16	128	128

注意黄的 Cb 和青的 Cr 都到了下限 16，蓝的 Cb 和红的 Cr 都到了上限 240。这正是第五节的结论：色差的极值出现在纯蓝、黄、纯红、青这几种颜色上。

一个直接的后果是：Y'CbCr 的整数组合有 $256^3$ 种，但其中只有一部分对应合法的 RGB 颜色。比如 $(Y, Cb, Cr) = (16, 240, 240)$ 又暗又蓝又红，换算回去 R、B 远大于 1、G 是负数。编码器在压缩时引入的误差也可能把颜色推出合法范围，所以解码后转 RGB 时必须截断到合法范围。

七、色度二次采样

1. 人眼对颜色细节不敏感

视觉研究早已发现，人眼对亮度变化的空间分辨能力，比对颜色变化的分辨能力高得多（Mullen，1985）：细密的黑白条纹能看清，同样细密的红绿条纹或蓝黄条纹，看起来就糊成一片了。

可以直接做实验（附录代码第 7 节）：把 Y'CbCr 的某些通道每 $k \times k$ 个像素平均一次，再插值回原尺寸，看看效果。

图 7：(b) 把两个色差通道做 8×8 的模糊，几乎看不出变化；(c) 只把亮度做同样的模糊，画面明显糊了。PSNR 是在整幅 R'G'B' 图像上算的。

模糊尺寸	只模糊色度：PSNR	剩下的数据量	只模糊亮度：PSNR	剩下的数据量
2×2	40.21 dB	50.0%	29.21 dB	75.0%
4×4	36.32 dB	37.5%	24.66 dB	68.8%
8×8	33.04 dB	34.4%	21.10 dB	67.2%

“剩下的数据量”是按样本数算的：只模糊色度时，Y' 全保留，两个色差各剩 $1/k^2$ ，总量是 $(1 + 2/k^2)/3$ 。同样丢掉一部分数据，模糊色度造成的误差小得多，而且看起来几乎没有区别。这就是色度二次采样（chroma subsampling）：降低色差通道的分辨率，亮度保持不变。

2. J:a:b 记法

二次采样的格式用三个数 J:a:b 表示，它描述一个宽 J 像素、高 2 行的区域：

J：区域宽度，通常是 4；
a：第一行里有几个色度样本；
b：第二行里新增几个色度样本（0 表示第二行和第一行共用）。

图 8：四种常见的色度二次采样格式。每个格子是一个像素，都有自己的亮度样本。

格式	色度分辨率	8 位时平均每像素	常见用途
4:4:4	不降	24 bit	专业制作、屏幕内容、需要精确颜色的场合
4:2:2	水平减半	16 bit	广播电视演播室、专业摄像机、ProRes
4:2:0	水平、垂直都减半	12 bit	几乎所有发行用的视频：网络视频、蓝光、电视广播
4:1:1	水平降到 1/4	12 bit	早期的 DV 磁带（NTSC）

4:2:0 的“0”很容易误解，它不是“第二行没有色度”，而是“第二行和第一行共用同一组色度”。所以 4:2:0 的每个色度样本对应 2×2 个像素，一幅 $W \times H$ 的图像，Cb 和 Cr 平面都是 $\frac{W}{2} \times \frac{H}{2}$ 。H.264 的 Main 和 High Profile 只支持 4:2:0（High 还允许只有亮度的黑白视频），4:2:2 和 4:4:4 要用专门的 High 4:2:2、High 4:4:4 Predictive Profile。4:2:0 也是本教程后面默认的格式。

3. 色度样本放在哪里

4:2:0 的一个色度样本代表 2×2 个亮度样本，但它的位置在哪里，不同标准的规定并不一样。

图 9：4:2:0 色度样本（橙色叉）相对于亮度样本（白圈）的三种位置。(b) 是 MPEG-2、H.264、H.265 的默认位置。

居中（图 9a）：JPEG 和 MPEG-1 的约定，色度在 4 个亮度样本正中间；
与左列对齐、上下居中（图 9b）：MPEG-2 起的视频标准默认采用，H.264、H.265 里叫 chroma_sample_loc_type = 0；
与左上角像素重合（图 9c）：chroma_sample_loc_type = 2。

H.264 和 H.265 可以在码流的 VUI 里标出 0～5 共 6 种位置（第 16 章）。位置搞错了，色度和亮度之间会错开不到半个像素，平时不易察觉，但在锐利的彩色边缘会出现色边。本章的代码为了简单，下采样用 2×2 平均、上采样用中心对齐的双线性插值，对应的是图 9(a) 的居中位置。

4. 下采样和上采样怎么做

4:2:0 的“损失”有多大，还取决于下采样和上采样用什么滤波器（附录代码第 8 节）：

下采样	上采样	R'G'B' PSNR	Cb 平面	Cr 平面
每 2×2 只取左上角一个	最近邻（复制）	36.57 dB	40.10 dB	40.68 dB
2×2 平均	最近邻（复制）	39.16 dB	43.02 dB	43.50 dB
2×2 平均	双线性插值	39.96 dB	43.68 dB	44.56 dB

和第 01 章的结论一样：下采样前要先低通滤波（平均就是最简单的低通），直接抽取会产生混叠，PSNR 低了 2.6 dB；上采样用插值比简单复制好。实际的转换工具（FFmpeg 的 swscale、libplacebo 等）用更长的滤波器，效果更好一些。

作为参照，只做 8 位量化、不做二次采样（4:4:4）时 PSNR 是 52.60 dB。4:2:0 把 PSNR 拉低到 40 dB 左右，但正如图 7 所示，这些误差在正常观看时几乎看不出来。

5. 色彩渗漏

二次采样也有看得出来的时候：颜色变化剧烈的细节。

图 10：(a)(b)(c) 宇航服上的任务徽章，4:2:0 之后文字边缘的颜色略微向外扩散，误差集中在彩色边缘；(d)(e) 宽 1、2、4 像素的红蓝条纹，1 像素宽的条纹在 4:2:0 之后完全变成了紫色。

1 像素宽的红蓝条纹，每个 2×2 块里一半红、一半蓝，色度一平均就成了紫色，原来的颜色信息彻底丢失。2 像素宽的条纹恰好与 2×2 块对齐，下采样没有损失，但双线性插值把相邻块的颜色混了进来，条纹变得发紫。4 像素宽的条纹只在交界处出现一列过渡色。

这类问题在自然视频里不明显，但在屏幕内容里很常见：红色的小字、彩色的细线、电子表格的网格线。所以录屏、远程桌面、游戏直播等场合常常要求 4:4:4，H.265 为此专门推出了“屏幕内容编码”扩展（第 24 章）。

八、常见的颜色错误

Y'CbCr 的数字本身不带“说明书”。同一组 $(Y, Cb, Cr)$ ，用不同的范围、不同的矩阵解释，会得到不同的颜色。编码时用一套、播放时用另一套，就会出事故。

图 11：同一份数据用不同方式解码。(b) 画面发灰，黑不黑、白不白；(c) 对比度过高，暗部和亮部大片截断；(d) 整体偏色，乍看不明显。

范围用错（附录代码第 9 节）：

有限范围当成全范围（图 11b）：黑（16）被当成 16/255 显示，白（235）被当成 235/255 显示。原图里的纯黑像素解码后平均是 16.0，纯白只到 235。画面像蒙了一层灰，PSNR 跌到 27.10 dB。这是最常见的视频颜色问题，网上“导出的视频发灰”大多是它；
全范围当成有限范围（图 11c）：码值被拉伸，暗部低于 16 的、亮部高于 235 的全部截断。被截成 0 或 255 的样本从原图的 11.2% 增加到 23.3%，PSNR 28.08 dB。

矩阵用错：BT.709 编码，按 BT.601 解码（图 11d），PSNR 35.44 dB。误差没有范围错误那么大，但很隐蔽：红、绿通道平均偏暗约 3 级，蓝通道偏亮 0.7 级，肤色微微发生变化。看纯色更明显：BT.709 的纯红按 BT.601 解码成 $(233,\ 0,\ 2)$ ，纯绿成了 $(20,\ 255,\ 9)$ 。高清视频被按标清的矩阵播放，或者反过来，是颜色“有点不对但说不上哪里不对”的常见原因。

怎么避免？在码流和文件里如实标记颜色参数，并在每一个环节检查它们（第十一节）。颜色参数没有标记时，各个软件只能自己猜，而且猜法不一样：FFmpeg 的 swscale 默认按 BT.601 换算，哪怕是高清视频；不少播放器（比如 mpv）则按分辨率猜，高清按 BT.709、标清按 BT.601。一边按 BT.601 生成、一边按 BT.709 播放，就是图 11(d) 的偏色。

九、误差预算与 YUV 文件

1. 每一步损失多少

把整条转换链串起来，看每一步的损失（附录代码第 10 节）：

图 12：从相机输出的 R'G'B' 到编码器的输入，经过的每一步和对应的 PSNR。

步骤	R'G'B' PSNR	与原图完全相同的像素
浮点 Y'CbCr 往返	无损（最大误差 $2.2\times10^{-16}$ ）	100%
8 位全范围 4:4:4	53.59 dB	33.6%
8 位有限范围 4:4:4	52.60 dB	26.0%
10 位有限范围 4:4:4	无损	100%
8 位有限范围 4:2:0	39.96 dB
10 位有限范围 4:2:0	40.18 dB

几个结论：

矩阵本身是可逆的，不损失信息；损失来自量化和二次采样；
8 位 Y'CbCr 存不下 8 位 R'G'B' 的全部信息，即使是 4:4:4，也只有约三成像素能原样还原。有限范围的码值更少，比全范围再低约 1 dB。10 位 4:4:4 则足以无损地还原 8 位 RGB；
在 4:2:0 面前，8 位和 10 位的差别几乎可以忽略（39.96 对 40.18 dB）：二次采样的损失比量化大得多。10 位视频的主要好处是减少色带和支持 HDR，而不是让 4:2:0 变得更准。

这些损失都发生在编码器开始工作之前。第 09 章讨论编码质量时，要记住参考图像本身已经是 4:2:0 了。

2. YUV 文件的布局

原始 YUV 文件和第 01 章的 RGB 文件一样，没有任何文件头，就是一串字节。同样是 4:2:0，字节的排列方式也有好几种：

图 13：一幅 4×2 的图像在 I420、NV12、YUY2 三种格式中的字节顺序。

I420（FFmpeg 叫 yuv420p）：平面格式。先存完整的 Y 平面，再存 U（Cb）平面，最后存 V（Cr）平面。软件编解码器（x264、libavcodec）最常用；
NV12：半平面格式。Y 平面之后，U、V 交错存在一个平面里。硬件编解码器、显卡、手机摄像头最常用；
YUY2（4:2:2）：打包格式。按 Y、U、Y、V 的顺序交错存放，每两个像素共用一对色度，摄像头的 USB 输出常用。

一帧 $W \times H$ 的 8 位 4:2:0 图像，大小是 $W \times H \times 1.5$ 字节。附录代码第 12 节把测试照片转成 I420，写出文件 astronaut_512x512_i420.yuv，大小是 393 216 字节，读回来逐字节核对无误。安装 FFmpeg 后（第 04 章），可以这样播放：

ffplay -f rawvideo -pixel_format yuv420p -video_size 512x512 -color_range tv -colorspace bt709 astronaut_512x512_i420.yuv

注意命令里专门告诉了播放器范围（tv 即有限范围）和矩阵（bt709），否则它只能去猜。

十、无损的颜色变换：YCoCg-R

Y'CbCr 的系数是无理数，8 位整数往返不可能无损（上一节已经看到）。有些场合需要严格无损：医学图像、屏幕内容、视频制作的中间文件。为此有一类只用整数加减和移位的颜色变换，最常用的是 YCoCg-R：

\begin{aligned} C_o &= R - B, & t &= B + \lfloor C_o / 2 \rfloor,\\ C_g &= G - t, & Y &= t + \lfloor C_g / 2 \rfloor \end{aligned}

逆变换把步骤倒过来：

t = Y - \lfloor C_g / 2 \rfloor,\quad G = C_g + t,\quad B = t - \lfloor C_o / 2 \rfloor,\quad R = B + C_o

每一步只是给某个变量加上（或减去）另一个变量的函数，倒过来做就能精确还原，不管那个函数里有没有取整。这种结构叫提升（lifting），第 07 章讲整数变换时还会用到。代价是 Co、Cg 的取值范围比原来大一倍，需要 9 位。

附录代码第 11 节在测试照片上验证：逆变换的最大误差是 0，完全无损。三个通道的一阶熵之和从 RGB 的 22.12 bit/像素降到 19.42 bit/像素，说明它同样起到了去相关的作用。H.265 的屏幕内容编码扩展和 H.266 都有类似的“自适应颜色变换”工具。

十一、码流里怎么标记颜色

为了避免第八节的事故，视频码流和容器里都有专门的字段说明颜色参数。H.264、H.265、H.266、AV1 共用 ITU-T H.273 定义的一套编号：

字段	含义	常见取值
`colour_primaries`	三原色和白点	1 = BT.709，5 = BT.601（625 行），6 = BT.601（525 行），9 = BT.2020
`transfer_characteristics`	传递函数	1 = BT.709，13 = sRGB，14 = BT.2020（10 位），16 = PQ，18 = HLG
`matrix_coefficients`	Y'CbCr 矩阵	0 = 不转换（直接存 GBR），1 = BT.709，5、6 = BT.601，9 = BT.2020
`video_full_range_flag`	范围	0 = 有限范围，1 = 全范围

在 H.264 里，这些字段放在 SPS 的 VUI 部分，第 16 章写码流解析器时会亲手把它们读出来。FFmpeg 对应的选项是 -color_primaries、-color_trc、-colorspace、-color_range。一个好习惯是：转码时显式写出这四个参数，不要依赖默认值。

十二、历史小插曲

年份	事件
1802	托马斯·杨（Thomas Young）提出人眼中有三种感光单元的假说，后由亥姆霍兹发展为三色理论
1861	麦克斯韦用红、绿、蓝三张滤色片拍摄并叠加投影，展示了世界上第一张彩色照片（一条格子花呢缎带）
1931	CIE 发布 1931 标准观察者和 XYZ 色彩空间
1953	NTSC 彩色电视：亮度加两个色差，与黑白电视兼容，亮度系数 0.299 / 0.587 / 0.114 由此而来
1976	CIE 发布 Lab* 和 u'v' 等更接近感知均匀的色彩空间
1982	CCIR 601（今 BT.601）规定数字视频的 Y'CbCr、4:2:2 采样和 16～235 的有限范围
1990	BT.709 首次发布，定义高清电视的原色和传递函数
1996	惠普和微软提出 sRGB，1999 年成为国际标准 IEC 61966-2-1
2012	BT.2020 发布，定义超高清电视的宽色域
2014	SMPTE ST 2084 定义 PQ 曲线
2016	BT.2100 发布，统一 HDR 的 PQ 和 HLG 两种方案

彩色电视当年面临的约束和今天的视频编码惊人地相似：频带有限，必须压缩；还要和已有的黑白电视兼容，所以必须单独传一个亮度信号。“亮度 + 色差 + 降低色差带宽”这个 1950 年代的方案，七十年后依然是每一个视频编码器的第一步。

小结

颜色是三维的：人眼只有三种视锥细胞，光谱被投影成三个数。CIE XYZ 是标准的三维坐标，其中 $Y$ 就是亮度。
色度图与色域： $(x,\ y)$ 去掉了亮度，所有颜色落在马蹄形内。设备的三原色构成色域三角形；BT.709 只覆盖 BT.2020 的约一半。
亮度系数由原色和白点决定：RGB → XYZ 矩阵可以从原色色度和白点推导，第二行就是 $K_R$ 、 $K_G$ 、 $K_B$ 。BT.709 是 0.2126 / 0.7152 / 0.0722，BT.601 的 0.299 / 0.587 / 0.114 来自 1953 年的 NTSC。
gamma 是为了人眼：人眼对亮度的感受近似立方根。8 位 gamma 编码的每一级 ΔL* 都不超过 0.51，线性编码要 11 位才能做到。HDR 用 PQ 或 HLG 代替传统 gamma。
Y'CbCr = 亮度 + 两个色差： $P_B = (B' - Y')/(2(1-K_B))$ ， $P_R = (R' - Y')/(2(1-K_R))$ 。它把约 89% 的能量集中到 Y'，与针对单张图像最优的 PCA 几乎一样好。
有限范围：8 位视频的黑是 16、白是 235，色差 16～240；JPEG 和屏幕录制常用全范围 0～255。
色度二次采样：人眼对颜色细节不敏感，4:2:0 把色度的分辨率减半，数据量从 24 降到 12 bit/像素，看起来几乎没有区别，但细小的彩色细节会渗色。
颜色参数必须标记：范围、矩阵、原色、传递函数四项，任何一项编码和解码不一致，都会造成发灰、过曝或偏色。

下一章开始讲压缩本身：四种冗余怎么分别消除，“预测 + 变换 + 量化 + 熵编码”的混合编码框架是怎么搭起来的，以及从 H.261 到 H.266 几十年来到底改进了什么。

练习题

（编码） 按 BT.601、8 位有限范围，把纯红 $(R', G', B') = (1, 0, 0)$ 编码成 $(Y, Cb, Cr)$ 。再按 BT.2020、10 位有限范围编码同一个纯红。
（推导） 证明 $P_B = (B' - Y')/(2(1-K_B))$ 的取值范围恰好是 $[-0.5,\ 0.5]$ ，并指出取到两个端点的颜色。
（码率） 1080p、60 帧/秒、10 位、4:2:2 的原始视频码率是多少？专业演播室常用的 3G-SDI 接口标称速率约 2.97 Gbit/s，装得下吗？
（gamma） 第四节说线性编码要 11 位才能让最大 ΔL* 不超过 8 位 sRGB 的 0.51。如果只要求每一级 ΔL* 不超过 1（刚好能察觉），线性编码需要几位？（提示：线性编码最大的台阶是第一级。）
（色域） 用第二节的表格，算出 DCI-P3 的色域面积是 BT.709 的多少倍（分别按 xy 和 u'v' 计算）。
（编程） 仿照附录代码第 8 节的 subsample_420，写一个 4:2:2（只在水平方向把色度减半）的版本，下采样用平均、上采样用双线性插值，计算 R'G'B' PSNR，并与 4:4:4 和 4:2:0 比较。
（YCoCg-R） 输入是 8 位 RGB 时， $C_o$ 和 $C_g$ 的理论取值范围各是多少？为什么需要 9 位？

参考答案

BT.601： $Y = 16 + 219 \times 0.299 = 81.5 \to 81$ ， $Cb = 128 + 224 \times (-0.1687) = 90.2 \to 90$ ， $Cr = 128 + 224 \times 0.5 = 240$ ，即 (81, 90, 240)。BT.2020、10 位： $Y = 4 \times (16 + 219 \times 0.2627) = 294.1 \to 294$ ， $Cb = 4 \times (128 - 224 \times 0.1396) = 386.9 \to 387$ ， $Cr = 4 \times 240 = 960$ ，即 (294, 387, 960)。两个结果都由附录代码的 rgb_to_ycbcr 实际算出。对比正文表格里 BT.709 的纯红 (63, 102, 240)，同一个“纯红”在不同标准下的 Y 值差别很大，这正是矩阵用错时颜色会变的原因。
$B' - Y' = (1-K_B)B' - K_R R' - K_G G'$ 。 $B'$ 的系数为正， $R'$ 、 $G'$ 的系数为负，且 $K_R + K_G = 1 - K_B$ 。取 $B' = 1$ 、 $R' = G' = 0$ （纯蓝）得最大值 $1 - K_B$ ；取 $B' = 0$ 、 $R' = G' = 1$ （黄）得最小值 $-(1 - K_B)$ 。除以 $2(1 - K_B)$ ，范围恰好是 $[-0.5,\ 0.5]$ 。同理， $P_R$ 在纯红处取 0.5，在青色处取 −0.5。
4:2:2 时每个像素平均 1 个亮度样本加 1 个色度样本（Cb、Cr 各占一半），10 位时 20 bit/像素。码率 $= 1920 \times 1080 \times 60 \times 20 = 2\,488\,320\,000$ bit/s，约 2.49 Gbit/s，小于 2.97 Gbit/s，装得下。实际上 3G-SDI 正是为 1080p60 4:2:2 10 位设计的，多出来的带宽用于行、场消隐期间的辅助数据（音频、时间码等）。
线性编码的码值正比于亮度，暗部 $L^* = 903.3\,Y$ ，所以第一级台阶 $\Delta L^* = 903.3 / (2^N - 1)$ 。要求它不超过 1，需要 $2^N - 1 \ge 903.3$ ： $N = 9$ 时为 $903.3/511 = 1.77$ ，不够； $N = 10$ 时为 $903.3/1023 = 0.88$ ，满足。即使标准放宽到“刚好能察觉”，线性编码也要 10 位才赶得上 8 位 gamma 编码。
xy： $0.1520 / 0.1120 \approx 1.36$ 倍；u'v'： $0.0815 / 0.0649 \approx 1.26$ 倍。在感知更均匀的 u'v' 坐标下，DCI-P3 比 BT.709 大约 26%。
参考实现：把每个色度平面 p 改成 np.round(p.reshape(512, 256, 2).mean(axis=2))（每行相邻两个样本平均），再用 resize_linear(small, p.shape) 插值回原尺寸。运行结果：4:2:2 的 R'G'B' PSNR 为 41.75 dB，介于 4:4:4 的 52.60 dB 和 4:2:0 的 39.96 dB 之间。只在一个方向减半，损失也大约只有 4:2:0 的一半。
$C_o = R - B$ ，取值从 $0 - 255 = -255$ 到 $255 - 0 = 255$ ，共 511 种可能。 $t = B + \lfloor C_o/2 \rfloor$ 介于 $B$ 和 $R$ 之间，所以在 0～255 内， $C_g = G - t$ 同样是 −255～255。511 种取值需要 9 位才能存下。 $Y = t + \lfloor C_g/2 \rfloor$ 介于 $t$ 和 $G$ 之间，仍在 0～255 内，8 位足够。附录代码在测试照片上实测的范围是 Co 从 −74 到 187、Cg 从 −112 到 64，都在理论范围内。

参考资料

C. Poynton. Digital Video and HD: Algorithms and Interfaces. 2nd ed. Morgan Kaufmann, 2012.（gamma、luma、Y'CbCr、色度二次采样最权威、最细致的参考书，本章大量观点出自此书）
E. Reinhard, E. A. Khan, A. O. Akyüz, G. M. Johnson. Color Imaging: Fundamentals and Applications. A K Peters, 2008.
G. Wyszecki, W. S. Stiles. Color Science: Concepts and Methods, Quantitative Data and Formulae. 2nd ed. Wiley, 1982.（色度学的经典大部头，CIE 配色函数的完整表格）
C. Wyman, P.-P. Sloan, P. Shirley. Simple Analytic Approximations to the CIE XYZ Color Matching Functions. Journal of Computer Graphics Techniques, 2(2): 1–11, 2013.
K. T. Mullen. The Contrast Sensitivity of Human Colour Vision to Red-Green and Blue-Yellow Chromatic Gratings. The Journal of Physiology, 359: 381–400, 1985.
H. S. Malvar, G. J. Sullivan. YCoCg-R: A Color Space with RGB Reversibility and Low Dynamic Range. Joint Video Team document JVT-I014, 2003.
标准文档：ITU-R BT.601、BT.709、BT.1886、BT.2020、BT.2087、BT.2100；SMPTE ST 2084；IEC 61966-2-1（sRGB）；ITU-T H.273（码流中的颜色参数编号）。

附录：完整代码

下面的代码复现了本章所有数值：配色函数与色度坐标、RGB → XYZ 矩阵与亮度系数、色域面积与色域转换矩阵、传递函数与 ΔL*、Y'CbCr 矩阵与彩条编码、去相关与主成分分析、色度与亮度模糊对比、色度二次采样、范围与矩阵用错、误差预算、YCoCg-R，以及 YUV 文件布局，最后画出只模糊色度与只模糊亮度的对比图。依赖 numpy、matplotlib 和 scikit-image，保存为 example.py 后运行 python example.py 即可。正文里的配图由同目录的 make_figures.py 生成。

"""视频编码自学教程 02：颜色与色彩空间 —— 附录代码

1. 配色函数与色度坐标：用解析近似公式算出 CIE 1931 配色函数和光谱轨迹；
2. 由三原色和白点推导 RGB -> XYZ 矩阵，亮度系数 0.2126 / 0.7152 / 0.0722 从哪里来；
3. 色域面积：BT.709、DCI-P3、BT.2020；BT.709 -> BT.2020 的转换矩阵；
4. 传递函数（gamma）：线性编码与 sRGB 编码每一级的感知步长 ΔL*，以及 PQ、HLG；
5. Y'CbCr：BT.601 / 709 / 2020 的矩阵，8 位有限范围的整数值（彩条）；
6. 去相关：RGB 与 Y'CbCr 三个通道的相关系数、方差占比，主成分分析；
7. 人眼对色度细节不敏感：只模糊色度与只模糊亮度的对比；
8. 色度二次采样：4:4:4 / 4:2:2 / 4:2:0 的数据量，4:2:0 的几种上下采样组合；
9. 常见错误：范围（有限 / 全范围）用错，矩阵（601 / 709）用错；
10. 误差预算：整条转换链每一步损失多少；
11. 无损的颜色变换 YCoCg-R；
12. YUV 文件的内存布局：I420、NV12、YUY2，写出并读回一个 I420 文件。
运行：python example.py（需要 numpy、matplotlib、scikit-image）
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data as skdata

plt.rcParams["font.sans-serif"] = ["Microsoft YaHei", "SimHei", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False


# ---------------- 1. 配色函数与色度坐标 ----------------
def _g(lam, mu, s1, s2):
    """左右宽度不同的高斯函数。"""
    s = np.where(lam < mu, s1, s2)
    return np.exp(-0.5 * ((lam - mu) / s) ** 2)


def cmf(lam):
    """CIE 1931 2° 标准观察者配色函数的解析近似（Wyman, Sloan, Shirley, 2013）。"""
    x = (1.056 * _g(lam, 599.8, 37.9, 31.0) + 0.362 * _g(lam, 442.0, 16.0, 26.7)
         - 0.065 * _g(lam, 501.1, 20.4, 26.2))
    y = 0.821 * _g(lam, 568.8, 46.9, 40.5) + 0.286 * _g(lam, 530.9, 16.3, 31.1)
    z = 1.217 * _g(lam, 437.0, 11.8, 36.0) + 0.681 * _g(lam, 459.0, 26.0, 13.8)
    return np.stack([x, y, z], axis=-1)


print("=== 1. 配色函数与色度坐标 ===")
lam = np.arange(380, 781, 1.0)
xyz_bar = cmf(lam)
for name, col in zip("xyz", xyz_bar.T):
    print(f"{name}bar 峰值在 {lam[np.argmax(col)]:.0f} nm，峰值 {col.max():.3f}")
for l0 in [450, 520, 580, 620]:
    X, Y, Z = cmf(np.array([float(l0)]))[0]
    print(f"{l0} nm 的单色光：色度坐标 x = {X / (X + Y + Z):.4f}，y = {Y / (X + Y + Z):.4f}")

# ---------------- 2. RGB -> XYZ 矩阵 ----------------
print("\n=== 2. 由三原色和白点推导 RGB -> XYZ 矩阵 ===")
PRIMARIES = {                                     # (R, G, B) 三原色的 xy 色度坐标
    "BT.709": [(0.640, 0.330), (0.300, 0.600), (0.150, 0.060)],
    "DCI-P3": [(0.680, 0.320), (0.265, 0.690), (0.150, 0.060)],
    "BT.2020": [(0.708, 0.292), (0.170, 0.797), (0.131, 0.046)],
    "NTSC 1953": [(0.67, 0.33), (0.21, 0.71), (0.14, 0.08)],
}
D65 = (0.3127, 0.3290)
ILL_C = (0.3101, 0.3162)                          # NTSC 1953 用的 C 光源白点


def xy_to_xyz(x, y):
    """亮度 Y = 1 时，色度坐标 (x, y) 对应的 XYZ。"""
    return np.array([x / y, 1.0, (1 - x - y) / y])


def rgb_to_xyz_matrix(prims, white):
    P = np.column_stack([xy_to_xyz(*p) for p in prims])   # 每列是一个原色（Y = 1）
    S = np.linalg.solve(P, xy_to_xyz(*white))             # 缩放，使 R = G = B = 1 时恰好是白点
    return P * S


np.set_printoptions(precision=4, suppress=True)
M709 = rgb_to_xyz_matrix(PRIMARIES["BT.709"], D65)
print("BT.709 RGB -> XYZ 矩阵：")
print(M709)
for name, white in [("BT.709", D65), ("BT.2020", D65), ("NTSC 1953", ILL_C)]:
    M = rgb_to_xyz_matrix(PRIMARIES[name], white)
    print(f"{name:<10s} 第二行（亮度系数 Kr, Kg, Kb）= {M[1]}")

# ---------------- 3. 色域 ----------------
print("\n=== 3. 色域面积 ===")


def tri_area(pts):
    (x1, y1), (x2, y2), (x3, y3) = pts
    return abs((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) / 2


def xy_to_uv(x, y):
    """CIE 1976 u'v' 色度坐标（比 xy 更接近感知均匀）。"""
    d = -2 * x + 12 * y + 3
    return 4 * x / d, 9 * y / d


for name in ["BT.709", "DCI-P3", "BT.2020"]:
    a_xy = tri_area(PRIMARIES[name])
    a_uv = tri_area([xy_to_uv(*p) for p in PRIMARIES[name]])
    r_xy = a_xy / tri_area(PRIMARIES["BT.2020"])
    r_uv = a_uv / tri_area([xy_to_uv(*p) for p in PRIMARIES["BT.2020"]])
    print(f"{name:<8s} xy 面积 {a_xy:.4f}（BT.2020 的 {r_xy:.1%}），u'v' 面积 {a_uv:.4f}（{r_uv:.1%}）")
M2020 = rgb_to_xyz_matrix(PRIMARIES["BT.2020"], D65)
M_709_to_2020 = np.linalg.solve(M2020, M709)
print("BT.709 RGB -> BT.2020 RGB（线性光）：")
print(M_709_to_2020)
print("BT.709 的纯红 (1, 0, 0) 在 BT.2020 里是", M_709_to_2020 @ [1, 0, 0])
print("BT.2020 的纯绿 (0, 1, 0) 在 BT.709 里是", np.linalg.solve(M_709_to_2020, [0, 1, 0]),
      "（有负值：超出 BT.709 色域）")

# ---------------- 4. 传递函数 ----------------
print("\n=== 4. 传递函数（gamma） ===")


def bt709_oetf(L):
    L = np.asarray(L, dtype=float)
    return np.where(L < 0.018, 4.5 * L, 1.099 * L ** 0.45 - 0.099)


def srgb_eotf(V):
    """sRGB 编码值 -> 线性光（0~1）。"""
    V = np.asarray(V, dtype=float)
    return np.where(V <= 0.04045, V / 12.92, ((V + 0.055) / 1.055) ** 2.4)


def lstar(Y):
    """CIE 1976 明度 L*（0~100），Y 为相对亮度（0~1）。"""
    Y = np.asarray(Y, dtype=float)
    return np.where(Y > (6 / 29) ** 3, 116 * np.cbrt(Y) - 16, Y * 903.2963)


def pq_inverse_eotf(nits):
    """SMPTE ST 2084（PQ）：绝对亮度（cd/m^2）-> 编码值 0~1。"""
    m1, m2 = 2610 / 16384, 2523 / 4096 * 128
    c1, c2, c3 = 3424 / 4096, 2413 / 4096 * 32, 2392 / 4096 * 32
    Y = np.asarray(nits, dtype=float) / 10000
    return ((c1 + c2 * Y ** m1) / (1 + c3 * Y ** m1)) ** m2


def hlg_oetf(E):
    a = 0.17883277
    b, c = 1 - 4 * a, 0.5 - a * np.log(4 * a)
    E = np.asarray(E, dtype=float)
    return np.where(E <= 1 / 12, np.sqrt(3 * E), a * np.log(np.maximum(12 * E - b, 1e-12)) + c)


print(f"18% 灰：明度 L* = {lstar(0.18):.1f}；BT.709 OETF 把它编码成 {bt709_oetf(0.18):.4f}，"
      f"L = 0.5 编码成 {bt709_oetf(0.5):.4f}")
codes = np.arange(256)
dL_lin = np.diff(lstar(codes / 255))                 # 线性编码：码值正比于亮度
dL_srgb = np.diff(lstar(srgb_eotf(codes / 255)))     # sRGB 编码
print(f"8 位线性编码：第一级 ΔL* = {dL_lin[0]:.2f}，最大 ΔL* = {dL_lin.max():.2f}，"
      f"ΔL* > 1 的台阶有 {np.sum(dL_lin > 1)} 个，"
      f"L* < 50 的暗部只分到 {np.sum(lstar(codes / 255) < 50)} 个码值")
print(f"8 位 sRGB 编码：第一级 ΔL* = {dL_srgb[0]:.2f}，最大 ΔL* = {dL_srgb.max():.2f}，"
      f"ΔL* > 1 的台阶有 {np.sum(dL_srgb > 1)} 个，"
      f"L* < 50 的暗部分到 {np.sum(lstar(srgb_eotf(codes / 255)) < 50)} 个码值")
for n in range(8, 17):
    c = np.arange(2 ** n)
    if np.diff(lstar(c / (2 ** n - 1))).max() <= dL_srgb.max():
        print(f"线性编码要 {n} 位，最大 ΔL* 才能不超过 8 位 sRGB 的 {dL_srgb.max():.2f}")
        break
for nits in [0.1, 1, 100, 1000, 10000]:
    print(f"PQ：{nits:>7g} cd/m^2 -> 编码值 {pq_inverse_eotf(nits):.4f}")
print(f"HLG：E = 1/12 -> {hlg_oetf(1 / 12):.4f}，E = 1 -> {hlg_oetf(1.0):.4f}")

# ---------------- 5. Y'CbCr ----------------
print("\n=== 5. Y'CbCr ===")
KRKB = {"BT.601": (0.299, 0.114), "BT.709": (0.2126, 0.0722), "BT.2020": (0.2627, 0.0593)}


def ycbcr_matrix(kr, kb):
    """R'G'B'（0~1）-> Y'（0~1）、Pb、Pr（-0.5~0.5）。"""
    kg = 1 - kr - kb
    return np.array([[kr, kg, kb],
                     [-kr / (2 * (1 - kb)), -kg / (2 * (1 - kb)), 0.5],
                     [0.5, -kg / (2 * (1 - kr)), -kb / (2 * (1 - kr))]])


for name, (kr, kb) in KRKB.items():
    print(f"{name}：Kr = {kr}，Kb = {kb}，Kg = {1 - kr - kb:.4f}")
    print(ycbcr_matrix(kr, kb))
A709 = ycbcr_matrix(*KRKB["BT.709"])
print("BT.709 逆矩阵（Y'PbPr -> R'G'B'）：")
print(np.linalg.inv(A709))


def rgb_to_ycbcr(rgb, std="BT.709", full=False, bits=8):
    """R'G'B'（0~1 浮点）-> 整数 Y'CbCr。"""
    ypp = rgb @ ycbcr_matrix(*KRKB[std]).T
    s = 2 ** (bits - 8)
    if full:
        scale, off = np.array([255, 255, 255]) * s, np.array([0, 128, 128]) * s
        hi = 2 ** bits - 1
    else:
        scale, off = np.array([219, 224, 224]) * s, np.array([16, 128, 128]) * s
        hi = 2 ** bits - 1
    return np.clip(np.round(ypp * scale + off), 0, hi).astype(np.int32)


def ycbcr_to_rgb(ycc, std="BT.709", full=False, bits=8):
    """整数 Y'CbCr -> R'G'B'（0~1 浮点，未截断）。"""
    s = 2 ** (bits - 8)
    if full:
        scale, off = np.array([255, 255, 255]) * s, np.array([0, 128, 128]) * s
    else:
        scale, off = np.array([219, 224, 224]) * s, np.array([16, 128, 128]) * s
    ypp = (ycc - off) / scale
    return ypp @ np.linalg.inv(ycbcr_matrix(*KRKB[std])).T


bars = {"白": (1, 1, 1), "黄": (1, 1, 0), "青": (0, 1, 1), "绿": (0, 1, 0),
        "品红": (1, 0, 1), "红": (1, 0, 0), "蓝": (0, 0, 1), "黑": (0, 0, 0)}
print("100% 彩条的 8 位 BT.709 有限范围编码 (Y, Cb, Cr)：")
for name, rgb in bars.items():
    print(f"  {name:<3s}{tuple(int(v) for v in rgb_to_ycbcr(np.array(rgb, float)))}")

# ---------------- 6. 去相关 ----------------
print("\n=== 6. 去相关：RGB 与 Y'CbCr ===")
rgb8 = skdata.astronaut()                            # 512x512x3，uint8
rgb = rgb8.astype(float) / 255
ypp = rgb @ A709.T                                   # Y'、Pb、Pr（浮点）
for name, arr, labels in [("R'G'B'", rgb, "RGB"), ("Y'CbCr", ypp, ["Y'", "Cb", "Cr"])]:
    flat = arr.reshape(-1, 3)
    C = np.corrcoef(flat.T)
    var = flat.var(axis=0)
    print(f"{name} 相关系数矩阵：")
    print(C)
    print(f"  方差占比：" + "，".join(f"{l} {v / var.sum():.1%}" for l, v in zip(labels, var)))
cov = np.cov(rgb.reshape(-1, 3).T)
w, V = np.linalg.eigh(cov)
w, V = w[::-1], V[:, ::-1]
pc1 = V[:, 0] * np.sign(V[:, 0].sum())
print(f"主成分分析：三个主成分的方差占比 {w[0] / w.sum():.1%}，{w[1] / w.sum():.1%}，{w[2] / w.sum():.1%}")
print(f"  第一主成分方向（归一化成和为 1）：{pc1 / pc1.sum()}，对比 BT.709 亮度系数 [0.2126 0.7152 0.0722]")

# ---------------- 7. 只模糊色度 vs 只模糊亮度 ----------------
print("\n=== 7. 只模糊色度 vs 只模糊亮度 ===")


def block_mean(img, k):
    h, w = img.shape
    return img.reshape(h // k, k, w // k, k).mean(axis=(1, 3))


def resize_linear(img, shape):
    """双线性缩放（像素中心对齐，边缘外推取边界值）。"""
    h, w = img.shape
    H, W = shape
    yi = np.clip((np.arange(H) + 0.5) * h / H - 0.5, 0, h - 1)
    xi = np.clip((np.arange(W) + 0.5) * w / W - 0.5, 0, w - 1)
    y0, x0 = np.floor(yi).astype(int), np.floor(xi).astype(int)
    y1, x1 = np.minimum(y0 + 1, h - 1), np.minimum(x0 + 1, w - 1)
    fy, fx = (yi - y0)[:, None], (xi - x0)[None, :]
    top = img[y0][:, x0] * (1 - fx) + img[y0][:, x1] * fx
    bot = img[y1][:, x0] * (1 - fx) + img[y1][:, x1] * fx
    return top * (1 - fy) + bot * fy


def blur(plane, k):
    return resize_linear(block_mean(plane, k), plane.shape)


def to8(rgb_float):
    return np.clip(np.round(rgb_float * 255), 0, 255)


def psnr(a, b, peak=255.0):
    mse = np.mean((np.asarray(a, float) - np.asarray(b, float)) ** 2)
    return 10 * np.log10(peak ** 2 / mse)


inv709 = np.linalg.inv(A709)
for k in [2, 4, 8]:
    chroma_blur = ypp.copy()
    chroma_blur[..., 1] = blur(ypp[..., 1], k)
    chroma_blur[..., 2] = blur(ypp[..., 2], k)
    luma_blur = ypp.copy()
    luma_blur[..., 0] = blur(ypp[..., 0], k)
    p_c = psnr(rgb8, to8(chroma_blur @ inv709.T))
    p_l = psnr(rgb8, to8(luma_blur @ inv709.T))
    print(f"{k}x{k} 模糊：只模糊色度 PSNR = {p_c:.2f} dB（数据量 {(1 + 2 / k ** 2) / 3:.1%}），"
          f"只模糊亮度 PSNR = {p_l:.2f} dB（数据量 {(1 / k ** 2 + 2) / 3:.1%}）")

# ---------------- 8. 色度二次采样 ----------------
print("\n=== 8. 色度二次采样 ===")
for name, a, b in [("4:4:4", 4, 4), ("4:2:2", 2, 2), ("4:2:0", 2, 0), ("4:1:1", 1, 1)]:
    chroma_per_4x2 = 2 * (a + b)                     # J = 4 的 4x2 区域里两路色度的样本数
    bpp = (8 * 8 + 8 * chroma_per_4x2) / 8           # 8 个亮度样本 + 色度样本，平均到每个像素
    print(f"{name}：每 4x2 个像素有 8 个 Y'、{chroma_per_4x2} 个色度样本，8 位时平均 {bpp:.0f} bit/像素")

ycc8 = rgb_to_ycbcr(rgb)                             # 8 位有限范围 4:4:4
ref = to8(ycbcr_to_rgb(ycc8))                        # 只经过 8 位量化的参考
print(f"4:4:4（只做 8 位有限范围量化）RGB PSNR = {psnr(rgb8, ref):.2f} dB")


def subsample_420(ycc, down="mean", up="linear"):
    out = ycc.astype(float).copy()
    for c in (1, 2):
        p = ycc[..., c].astype(float)
        small = block_mean(p, 2) if down == "mean" else p[::2, ::2]
        small = np.round(small)                      # 色度平面同样存成 8 位整数
        if up == "linear":
            out[..., c] = resize_linear(small, p.shape)
        else:
            out[..., c] = np.repeat(np.repeat(small, 2, axis=0), 2, axis=1)
    return out


for down, up in [("drop", "nearest"), ("mean", "nearest"), ("mean", "linear")]:
    rec = subsample_420(ycc8, down, up)
    p_rgb = psnr(rgb8, to8(ycbcr_to_rgb(np.round(rec))))
    p_cb = psnr(ycc8[..., 1], np.round(rec[..., 1]))
    p_cr = psnr(ycc8[..., 2], np.round(rec[..., 2]))
    print(f"4:2:0 下采样 {down:<5s} 上采样 {up:<8s}：RGB PSNR = {p_rgb:.2f} dB，"
          f"Cb 平面 {p_cb:.2f} dB，Cr 平面 {p_cr:.2f} dB")

# ---------------- 9. 常见错误 ----------------
print("\n=== 9. 范围与矩阵用错 ===")
lim = rgb_to_ycbcr(rgb)                              # 正确：BT.709 有限范围
wrong_full = to8(ycbcr_to_rgb(lim, full=True))       # 有限范围被当成全范围解码
print(f"有限范围当成全范围解码：PSNR = {psnr(rgb8, wrong_full):.2f} dB，"
      f"原图纯黑像素解码成 {wrong_full[rgb8.sum(axis=2) == 0].mean():.1f}（应为 0），"
      f"纯白（255,255,255）解码成 {ycbcr_to_rgb(np.array([235, 128, 128]), full=True)[0] * 255:.1f}")
full = rgb_to_ycbcr(rgb, full=True)
wrong_lim = to8(ycbcr_to_rgb(full, full=False))      # 全范围被当成有限范围解码
clipped = np.mean((wrong_lim == 0) | (wrong_lim == 255))
print(f"全范围当成有限范围解码：PSNR = {psnr(rgb8, wrong_lim):.2f} dB，"
      f"被截到 0 或 255 的样本占 {clipped:.1%}（原图 {np.mean((rgb8 == 0) | (rgb8 == 255)):.1%}）")
wrong_mat = to8(ycbcr_to_rgb(lim, std="BT.601"))     # BT.709 编码，BT.601 解码
print(f"BT.709 编码、BT.601 解码：PSNR = {psnr(rgb8, wrong_mat):.2f} dB，"
      f"各通道平均偏差 R {np.mean(wrong_mat[..., 0] - rgb8[..., 0]):+.2f}，"
      f"G {np.mean(wrong_mat[..., 1] - rgb8[..., 1]):+.2f}，B {np.mean(wrong_mat[..., 2] - rgb8[..., 2]):+.2f}")
for name in ["红", "绿", "蓝"]:
    c = rgb_to_ycbcr(np.array(bars[name], float))
    print(f"  纯{name} {bars[name]} 按 BT.601 解码成 {tuple(int(v) for v in to8(ycbcr_to_rgb(c, std='BT.601')))}")

# ---------------- 10. 误差预算 ----------------
print("\n=== 10. 误差预算 ===")
roundtrip = (rgb @ A709.T) @ inv709.T
print(f"浮点 Y'CbCr 往返：最大误差 {np.abs(roundtrip - rgb).max():.1e}")
for label, kw in [("8 位全范围 4:4:4", dict(full=True)), ("8 位有限范围 4:4:4", dict()),
                  ("10 位有限范围 4:4:4", dict(bits=10))]:
    rec = to8(ycbcr_to_rgb(rgb_to_ycbcr(rgb, **kw), **kw))
    same = np.mean(np.all(rec == rgb8, axis=2))
    p = "无穷大" if same == 1 else f"{psnr(rgb8, rec):.2f} dB"
    print(f"{label:<16s} RGB PSNR = {p}，与原图完全相同的像素占 {same:.1%}")
for label, bits in [("8 位有限范围 4:2:0", 8), ("10 位有限范围 4:2:0", 10)]:
    ycc = rgb_to_ycbcr(rgb, bits=bits)
    rec = to8(ycbcr_to_rgb(np.round(subsample_420(ycc)), bits=bits))
    print(f"{label:<16s} RGB PSNR = {psnr(rgb8, rec):.2f} dB")

# ---------------- 11. YCoCg-R ----------------
print("\n=== 11. 无损颜色变换 YCoCg-R ===")
R, G, B = (rgb8[..., i].astype(np.int32) for i in range(3))
Co = R - B
t = B + (Co >> 1)
Cg = G - t
Y = t + (Cg >> 1)
t2 = Y - (Cg >> 1)
G2 = Cg + t2
B2 = t2 - (Co >> 1)
R2 = B2 + Co
print(f"逆变换后最大误差：{max(np.abs(R2 - R).max(), np.abs(G2 - G).max(), np.abs(B2 - B).max())}"
      f"（完全无损）")
print(f"取值范围：Y {Y.min()}~{Y.max()}，Co {Co.min()}~{Co.max()}，Cg {Cg.min()}~{Cg.max()}"
      "（Co、Cg 需要 9 位）")


def entropy(v):
    _, counts = np.unique(v, return_counts=True)
    p = counts / counts.sum()
    return -np.sum(p * np.log2(p))


h_rgb = sum(entropy(c) for c in (R, G, B))
h_ycocg = sum(entropy(c) for c in (Y, Co, Cg))
print(f"三个通道的一阶熵之和：RGB {h_rgb:.2f} bit/像素，YCoCg-R {h_ycocg:.2f} bit/像素")

# ---------------- 12. YUV 文件布局 ----------------
print("\n=== 12. YUV 文件的内存布局 ===")
w, h = 4, 2
Yl = np.array([f"Y{i}" for i in range(w * h)]).reshape(h, w)
Ul = np.array([f"U{i}" for i in range(w // 2 * h // 2)]).reshape(h // 2, w // 2)
Vl = np.array([f"V{i}" for i in range(w // 2 * h // 2)]).reshape(h // 2, w // 2)
i420 = list(Yl.ravel()) + list(Ul.ravel()) + list(Vl.ravel())
nv12 = list(Yl.ravel()) + [s for uv in zip(Ul.ravel(), Vl.ravel()) for s in uv]
U2 = np.array([f"U{i}" for i in range(w // 2 * h)]).reshape(h, w // 2)    # 4:2:2：每行都有色度
V2 = np.array([f"V{i}" for i in range(w // 2 * h)]).reshape(h, w // 2)
yuy2 = [s for r in range(h) for c in range(w // 2)
        for s in (Yl[r, 2 * c], U2[r, c], Yl[r, 2 * c + 1], V2[r, c])]
print("4x2 的图像：")
print("  I420（YUV420p）：", " ".join(i420))
print("  NV12           ：", " ".join(nv12))
print("  YUY2（4:2:2）  ：", " ".join(yuy2))

ycc = rgb_to_ycbcr(rgb)
Yp = ycc[..., 0].astype(np.uint8)
Up = np.round(block_mean(ycc[..., 1].astype(float), 2)).astype(np.uint8)
Vp = np.round(block_mean(ycc[..., 2].astype(float), 2)).astype(np.uint8)
fname = "astronaut_512x512_i420.yuv"
with open(fname, "wb") as f:
    for plane in (Yp, Up, Vp):
        f.write(plane.tobytes())
raw = np.fromfile(fname, dtype=np.uint8)
W, H = 512, 512
Yr = raw[:W * H].reshape(H, W)
Ur = raw[W * H: W * H * 5 // 4].reshape(H // 2, W // 2)
Vr = raw[W * H * 5 // 4:].reshape(H // 2, W // 2)
print(f"写出 {fname}：{raw.size:,} 字节 = 512×512×1.5；读回后三个平面完全一致："
      f"{np.array_equal(Yr, Yp) and np.array_equal(Ur, Up) and np.array_equal(Vr, Vp)}")
print("播放命令：ffplay -f rawvideo -pixel_format yuv420p -video_size 512x512 "
      f"-color_range tv -colorspace bt709 {fname}")

# ---------------- 画图：只模糊色度 vs 只模糊亮度 ----------------
k = 8
cb = ypp.copy()
cb[..., 1], cb[..., 2] = blur(ypp[..., 1], k), blur(ypp[..., 2], k)
lb = ypp.copy()
lb[..., 0] = blur(ypp[..., 0], k)
fig, axes = plt.subplots(1, 3, figsize=(12, 4.4))
for ax, img, title in [(axes[0], rgb8, "原图"),
                       (axes[1], to8(cb @ inv709.T) / 255, "只模糊色度（8x8）"),
                       (axes[2], to8(lb @ inv709.T) / 255, "只模糊亮度（8x8）")]:
    ax.imshow(img)
    ax.set_title(title)
    ax.axis("off")
fig.tight_layout()
plt.show()

运行后的关键输出（节选）：

=== 2. 由三原色和白点推导 RGB -> XYZ 矩阵 ===
BT.709     第二行（亮度系数 Kr, Kg, Kb）= [0.2126 0.7152 0.0722]
BT.2020    第二行（亮度系数 Kr, Kg, Kb）= [0.2627 0.678  0.0593]
NTSC 1953  第二行（亮度系数 Kr, Kg, Kb）= [0.2989 0.5866 0.1144]
=== 3. 色域面积 ===
BT.709   xy 面积 0.1120（BT.2020 的 52.9%），u'v' 面积 0.0649（58.0%）
DCI-P3   xy 面积 0.1520（BT.2020 的 71.7%），u'v' 面积 0.0815（72.9%）
BT.2020 的纯绿 (0, 1, 0) 在 BT.709 里是 [-0.5876  1.1329 -0.1006] （有负值：超出 BT.709 色域）
=== 4. 传递函数（gamma） ===
18% 灰：明度 L* = 49.5；BT.709 OETF 把它编码成 0.4090，L = 0.5 编码成 0.7055
8 位线性编码：第一级 ΔL* = 3.54，最大 ΔL* = 3.54，ΔL* > 1 的台阶有 15 个，L* < 50 的暗部只分到 47 个码值
8 位 sRGB 编码：第一级 ΔL* = 0.27，最大 ΔL* = 0.51，ΔL* > 1 的台阶有 0 个，L* < 50 的暗部分到 119 个码值
线性编码要 11 位，最大 ΔL* 才能不超过 8 位 sRGB 的 0.51
PQ：    100 cd/m^2 -> 编码值 0.5081
PQ：   1000 cd/m^2 -> 编码值 0.7518
=== 6. 去相关：RGB 与 Y'CbCr ===
  方差占比：R 36.1%，G 31.5%，B 32.5%
  方差占比：Y' 88.6%，Cb 3.1%，Cr 8.3%
主成分分析：三个主成分的方差占比 89.3%，10.0%，0.7%
=== 7. 只模糊色度 vs 只模糊亮度 ===
8x8 模糊：只模糊色度 PSNR = 33.04 dB（数据量 34.4%），只模糊亮度 PSNR = 21.10 dB（数据量 67.2%）
=== 8. 色度二次采样 ===
4:4:4（只做 8 位有限范围量化）RGB PSNR = 52.60 dB
4:2:0 下采样 mean  上采样 linear  ：RGB PSNR = 39.96 dB，Cb 平面 43.68 dB，Cr 平面 44.56 dB
=== 9. 范围与矩阵用错 ===
有限范围当成全范围解码：PSNR = 27.10 dB，原图纯黑像素解码成 16.0（应为 0），纯白（255,255,255）解码成 235.0
全范围当成有限范围解码：PSNR = 28.08 dB，被截到 0 或 255 的样本占 23.3%（原图 11.2%）
BT.709 编码、BT.601 解码：PSNR = 35.44 dB，各通道平均偏差 R -3.17，G -3.10，B +0.72
=== 10. 误差预算 ===
8 位有限范围 4:4:4    RGB PSNR = 52.60 dB，与原图完全相同的像素占 26.0%
10 位有限范围 4:4:4   RGB PSNR = 无穷大，与原图完全相同的像素占 100.0%
8 位有限范围 4:2:0    RGB PSNR = 39.96 dB
10 位有限范围 4:2:0   RGB PSNR = 40.18 dB
=== 11. 无损颜色变换 YCoCg-R ===
逆变换后最大误差：0（完全无损）
三个通道的一阶熵之和：RGB 22.12 bit/像素，YCoCg-R 19.42 bit/像素
=== 12. YUV 文件的内存布局 ===
写出 astronaut_512x512_i420.yuv：393,216 字节 = 512×512×1.5；读回后三个平面完全一致：True