11.2 抖动与蓝噪声
站内 /graphics-anti-aliasing-methods/ 那篇在时间域上做抗锯齿——多采几帧再平均。本节换到空间域，而且一分采样都不加：输入还是那张 128×128 的 0..255 水平斜坡，输出还是 8 位，只是量化到 8 级。问题出在量化本身——步长 36，于是斜坡被切成一条条 36 级的色带，这就是带状伪影（banding）。
抖动（dithering）的全部作用是：把量化误差从低频搬到高频。它不减少误差的总量（误差预算是量化步长决定的，抖动改不了），改的是误差的频谱分布——人眼对低频的带状轮廓极其敏感，对高频的细密噪点则宽容得多。所以「抖动好不好」不能用 MAE 一个数判断，必须看误差落在哪个频段。
四个方案，一条斜坡
| 方案 | 阈值来源 | 误差落在哪 |
|---|---|---|
none | 常数 128（四舍五入） | 每条带内线性变化，低频、可见 |
ordered | Bayer 4×4（整数递推构造） | 高频，但图案规则，有谱峰 |
blue | void-and-cluster 16×16 rank 矩阵 | 高频且各向同性，低频功率最小 |
fs | 无阈值矩阵，误差向后扩散 | 高频、不规则，逐像素串行 |
Bayer 矩阵的整数递推是本节第一个可断言的东西：M_2n = [[4M+0, 4M+2], [4M+3, 4M+1]]，从 [[0]] 出发翻倍四次就是 16×16，元素恰好是 0..n²-1 的一个排列。量化则把阈值折进偏移量：off = (thr * step) // 256，q = clip((v + off) // step, 0, levels-1) * step。误差扩散走另一条路——它不查表，把每个像素的量化误差按 7/16、3/16、5/16、1/16 的权重分给右、左下、下、右下四个邻居。
S 层的三段核心（完整脚本见工作区 sim/cg_shade_11_2_dither.py）：
# snippet: 11.2 的三种抖动构造（节选）：Bayer 整数递推 / 带阈值量化 / Floyd-Steinberg 整数扩散
def bayer(n):
"""Bayer 矩阵的整数递推：M_2n = [[4M+0, 4M+2], [4M+3, 4M+1]]。"""
m = np.zeros((1, 1), dtype=np.int64)
k = 1
while k < n:
m = np.block([[4 * m + 0, 4 * m + 2], [4 * m + 3, 4 * m + 1]])
k *= 2
return m
def quantize(v, levels, thr):
"""带阈值的均匀量化：`thr` 为 Q8 阈值（0..255），返回量化后的 0..255 整数。"""
step = 255 // (levels - 1)
off = (thr * step) // 256
qi = np.clip((v + off) // step, 0, levels - 1)
return qi * step
def fs_dither(v, levels):
"""Floyd-Steinberg 误差扩散（整数）。误差缓冲以 16 为分母，避免浮点。"""
h, w = v.shape
step = 255 // (levels - 1)
out = np.zeros_like(v)
err = np.zeros((h + 1, w + 2), dtype=np.int64)
for y in range(h):
for x in range(w):
p = v[y, x] + err[y, x + 1] // 16
qi = min(levels - 1, max(0, (p + step // 2) // step))
q = qi * step
out[y, x] = q
e = p - q
err[y, x + 2] += e * 7
err[y + 1, x] += e * 3
err[y + 1, x + 1] += e * 5
err[y + 1, x + 2] += e
return out
err[y, x+1] // 16 是这里唯一会碰到负数的一行——误差有正有负，Python 的 // 向 -∞ 取整。搬到 TypeScript 时这一行必须写成 Math.floor，写成 | 0 或 Math.trunc 就会在负误差上差 1，累积几百个像素后整幅图都变了。这是本章第二次踩同一个坑（11.1 是 >> 的 32 位域）。
实测（/opt/homebrew/bin/python3 sim/cg_shade_11_2_dither.py）：
画布 / 级数 : 128x128 量化到 8 级（步长 36）
方案 MAE 低频(带状) 高频(噪声) 谱峰 有符号均值
none 8.875 6.4844 2.4531 91120.04 0.9688
ordered 11.8281 1.2922 11.2155 145152.0 -0.0156
blue 11.7617 1.3262 11.9453 32487.81 -0.1035
fs 11.4558 1.1705 11.2203 22720.98 -0.4507
低频功率占比 : blue 16x16 = 0.000217 white 16x16 = 0.052722 （蓝噪声把低频功率压低 243 倍）
这张表要横着读。总误差那一列（MAE）反而变大了：无抖动 8.875，三种抖动都在 11.5～11.8。这不是 bug——抖动的阈值偏移让误差不再逐像素抵消，MAE 上升是必然的代价。真正变化的是第二列：低频（带状）分量从 6.4844 掉到 1.2 上下，降了五倍；对应的第三列高频（噪声）分量从 2.4531 涨到 11～12。误差被搬走了，没被消灭。
有符号均值一列是「直流偏置」的体检：无抖动时是 0.9688（带内误差系统性偏正，这正是带状轮廓的来源），三种抖动都压到 0.1 以内。抖动的阈值矩阵必须均值为 128，否则它会引入一条整体色偏——这一条在 P 层被写成了断言。
最后一行是蓝噪声的定义性指标，也是本节唯一一个「谁更好」可以复算的地方：把 16×16 蓝噪声 rank 矩阵与同尺寸白噪声矩阵都做零均值化、求二维功率谱，再看半径小于 3 的低频段占总功率的比例。蓝噪声 0.000217，白噪声 0.052722——低两个数量级，243 倍。这正是「蓝」的含义：能量被推到高频去了。
P 层：Bayer、误差扩散与朴素 DFT
P 层用 TypeScript 重写同一套东西（core/src/lantern/cg_shade_11_dither.ts），要求四个量化平面与 Python 逐字节相同：
// file: lantern/core/src/lantern/cg_shade_11_dither.ts
// tier: P
// verify: node --test 'core/test/cg_shade_11_2_dither.test.ts'
// 11.2 抖动与蓝噪声的 P 层：Bayer 整数递推、带阈值量化、Floyd–Steinberg 整数误差扩散，
// 以及低频功率占比（蓝噪声的定义性指标）。与 Python 侧 sim/cg_shade_11_2_dither.py
// 做**同一件事**，两侧是两份独立代码，结果必须逐字节相同。
//
// 三个「Python 与 JS 语义不同」的坑，本节全部踩在明处：
// ① 误差扩散里 err 可为负，Python `//` 向 -inf 取整，故 JS 侧必须用 `Math.floor`，不能用 `| 0`；
// ② 量化偏移 `(thr*step)//256` 的分子非负，floor 与截断同义；
// ③ 谱分析用朴素 DFT（16x16 只有 65536 次乘加），避免引入 FFT 依赖。
/** Bayer 矩阵的整数递推：M_2n = [[4M+0, 4M+2], [4M+3, 4M+1]]。n 必须是 2 的幂。 */
export function bayer(n: number): number[][] {
let m: number[][] = [[0]];
let k = 1;
while (k < n) {
const size = k * 2;
const next: number[][] = Array.from({ length: size }, () => new Array<number>(size).fill(0));
for (let y = 0; y < k; y += 1) {
for (let x = 0; x < k; x += 1) {
const v = m[y][x] * 4;
next[y][x] = v;
next[y][x + k] = v + 2;
next[y + k][x] = v + 3;
next[y + k][x + k] = v + 1;
}
}
m = next;
k = size;
}
return m;
}
/**
* 把 n x n 阈值矩阵平铺到 w x h。
* `mode = "bayer"`：阈值 = r*16 + 8（r 是 0..n²-1 的序号，把阈值均分到 0..255）；
* `mode = "rank"`：阈值 = r（蓝噪声 rank 矩阵本身就是 0..255 的阈值）。
*/
export function tileThreshold(mat: number[][], w: number, h: number, mode: "bayer" | "rank"): Int32Array {
const n = mat.length;
const out = new Int32Array(w * h);
for (let y = 0; y < h; y += 1) {
for (let x = 0; x < w; x += 1) {
const r = mat[y % n][x % n];
out[y * w + x] = mode === "bayer" ? r * 16 + 8 : r;
}
}
return out;
}
/** 带阈值的均匀量化：q = clip((v + (thr*step)//256) // step, 0, levels-1) * step。 */
export function quantize(v: Int32Array, levels: number, thr: Int32Array): Int32Array {
const step = Math.floor(255 / (levels - 1));
const out = new Int32Array(v.length);
for (let i = 0; i < v.length; i += 1) {
const off = Math.floor((thr[i] * step) / 256);
let qi = Math.floor((v[i] + off) / step);
if (qi < 0) qi = 0;
else if (qi > levels - 1) qi = levels - 1;
out[i] = qi * step;
}
return out;
}
/** Floyd–Steinberg 误差扩散（整数）。误差缓冲以 16 为分母，全程无浮点。 */
export function fsDither(v: Int32Array, w: number, h: number, levels: number): Int32Array {
const step = Math.floor(255 / (levels - 1));
const out = new Int32Array(w * h);
const EW = w + 2;
const err = new Int32Array((h + 1) * EW);
for (let y = 0; y < h; y += 1) {
for (let x = 0; x < w; x += 1) {
// err 可为负 -> 必须 Math.floor（Python 的 // 向 -inf 取整）
const p = v[y * w + x] + Math.floor(err[y * EW + x + 1] / 16);
const qi = Math.min(levels - 1, Math.max(0, Math.floor((p + (step >> 1)) / step)));
const q = qi * step;
out[y * w + x] = q;
const e = p - q;
err[y * EW + x + 2] += e * 7;
err[(y + 1) * EW + x] += e * 3;
err[(y + 1) * EW + x + 1] += e * 5;
err[(y + 1) * EW + x + 2] += e;
}
}
return out;
}
/** 5x5 盒式低通（整数，边缘用 reflect 填充）。 */
export function box5(a: Int32Array, w: number, h: number): Int32Array {
const at = (x: number, y: number) => {
const xx = x < 0 ? -x : (x >= w ? 2 * w - 2 - x : x);
const yy = y < 0 ? -y : (y >= h ? 2 * h - 2 - y : y);
return a[yy * w + xx];
};
const out = new Int32Array(w * h);
for (let y = 0; y < h; y += 1) {
for (let x = 0; x < w; x += 1) {
let s = 0;
for (let dy = -2; dy <= 2; dy += 1) for (let dx = -2; dx <= 2; dx += 1) s += at(x + dx, y + dy);
out[y * w + x] = Math.floor(s / 25);
}
}
return out;
}
/** 低频（带状）误差：|5x5 盒式低通(量化误差)| 的均值。抖动要做的事就是把它压低。 */
export function bandingLf(q: Int32Array, v: Int32Array, w: number, h: number): number {
const e = new Int32Array(q.length);
for (let i = 0; i < q.length; i += 1) e[i] = q[i] - v[i];
const lf = box5(e, w, h);
let s = 0;
for (let i = 0; i < lf.length; i += 1) s += Math.abs(lf[i]);
return s / lf.length;
}
/** 二维功率谱（朴素 DFT，O(n^4)）。仅用于 16x16 的阈值矩阵。 */
export function powerSpectrum(mat: number[][]): number[][] {
const n = mat.length;
const out: number[][] = [];
for (let u = 0; u < n; u += 1) {
out.push(new Array<number>(n).fill(0));
for (let v = 0; v < n; v += 1) {
let re = 0;
let im = 0;
for (let y = 0; y < n; y += 1) {
for (let x = 0; x < n; x += 1) {
const ang = (-2 * Math.PI * (u * y + v * x)) / n;
re += mat[y][x] * Math.cos(ang);
im += mat[y][x] * Math.sin(ang);
}
}
out[u][v] = re * re + im * im;
}
}
return out;
}
/** 低频功率占比：半径 < cut 的频段功率 / 全频段功率。蓝噪声的定义性指标，越小越「蓝」。 */
export function lfPowerRatio(mat: number[][], cut = 3): number {
const n = mat.length;
let mean = 0;
for (const row of mat) for (const v of row) mean += v;
mean /= n * n;
const centered = mat.map((row) => row.map((v) => v - mean));
const P = powerSpectrum(centered);
let total = 0;
let low = 0;
for (let u = 0; u < n; u += 1) {
for (let v = 0; v < n; v += 1) {
const p = P[u][v];
total += p;
const ru = Math.min(u, n - u);
const rv = Math.min(v, n - v);
if (Math.hypot(ru, rv) < cut) low += p;
}
}
return low / total;
}
谱分析那一段刻意用朴素 DFT而不是 FFT：16×16 的矩阵只有 256 个频点、每个频点 256 次乘加，一共 65536 次——写成一团四重循环也就几毫秒，换来的是「不引入任何依赖」。蓝噪声的判据不需要快，需要的是能被独立复算。
node --test 的结果：
✔ Bayer 整数递推：M_2 与 M_4 逐项正确
✔ 四个量化平面与 Python 侧逐字节相同
✔ 低频（带状）误差与 S 层一致：无抖动最高，三种抖动都被压低
✔ 蓝噪声的定义性指标：低频功率占比比同尺寸白噪声低两个数量级
✔ 负例：把蓝噪声阈值换成常数 128 -> 退化为无抖动，低频必须回升
tests 5 / pass 5 / fail 0
第二条是跨语言断言：四个平面各 16384 字节，Python 与 TypeScript 逐字节相同——只要 Math.floor 那一行写错，fs 平面立刻对不上。第五条是负例：把蓝噪声的 rank 阈值换成常数 128，图案退化成无抖动，低频从 1.33 涨回 6.48，判据必须拒绝它。
小结
抖动是一类「用误差的空间分布换主观质量」的技术，它的判据必须分频段看：总误差（MAE）可以变大，低频（带状）必须变小，高频（噪声）变大是代价而不是缺陷。有序抖动与误差扩散都能把低频压到 1/5，区别在图案的规则度——Bayer 的谱峰是 145152，误差扩散只有 22720，肉眼看到的分别是「规则网格」和「随机颗粒」。蓝噪声介于两者之间：图案不规则（谱峰 32487），但低频功率最低（0.000217）。下一节把「谁更好」写成一条可复算的判据，而不是形容词。
阅读导航：上一节：11.1 时间抗锯齿 · 下一节：11.3 抗锯齿的度量 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。