字半色调原理详解（四）

12.3 图案抖动

12.4 非周期抖动

12.4.1 弗洛伊德‑斯坦伯格算法

彩色图像的抖动处理

12.4.2 基于空间填充曲线的抖动算法

像素遍历排序

基于皮亚诺曲线的图像区域划分

自适应分割

自适应分割计算流程

抖动单元量化

黑点簇位置微调

完整抖动算法流程

12.3 图案抖动

有序抖动常常会和一种名为图案抖动（Pattern Dithering）的技术相混淆。该方法将图像划分为若干大小为 N×N 的抖动单元，这一点和有序抖动的处理方式类似。每个抖动单元一共可以生成 (N^2+1) 个量化灰度等级。

算法执行流程如下：使用这 (N^2+1) 个量化等级作为索引，读取一张抖动图案查找表；针对每一个单元，利用式 (12.1) 计算单元内图像的平均亮度，再将该平均值量化至 (N^2+1) 个等级中的某一级。最后以量化后的等级值为索引，取出查找表中对应的抖动图案。

需要注意：如果抖动图案查找表直接选用 12.2 节算出的有序抖动所对应的 (N^2+1) 套图案，那么图案抖动算法生成的效果与有序抖动相近但并不完全等同。事实上，对于同一个平均亮度，图案抖动永远查表输出完全相同的图案；而有序抖动则会根据抖动单元内部像素亮度的实际分布，生成不一样的二值点阵图案。

抖动图案查找表的选取自由度很高。正因如此，图案抖动可用来为图像生成多种特殊渲染效果。该技术一个经典的实例就是行式打印机输出的图像（见图 12‑18）：这类效果的实现原理就是构造一张抖动表，让每一个量化等级都对应一组精心挑选、可重叠打印的字符。

即便在这类特殊效果场景下，基于阈值函数的抖动（有序抖动）仍然比图案抖动灵活得多。再次提醒读者：有序抖动生成的点簇形状由阈值函数的等值线几何形态决定，而阈值函数本身可以设计成任意形状。

12.4 非周期抖动

在各类非周期抖动算法当中，最早出现的一批算法包含随机调制抖动；不过现如今，该算法已经大多被效果更优的方案所取代。而另一种问世稍晚、现已成为经典且广为流行的算法便是弗洛伊德‑斯坦伯格（Floyd‑Steinberg）误差扩散算法。该算法最初面向灰度图像提出，也能够很方便地拓展应用于彩色图像。

12.4.1 弗洛伊德‑斯坦伯格算法

弗洛伊德‑斯坦伯格算法会计算当前像素量化后产生的误差，并将该误差补偿分摊到周边相邻像素上，以此将整幅图像的总量化误差降到最低。具体来说，假设我们要对由图像函数 (f(x,y)) 定义的图像执行 1‑bit 二值量化。从坐标 ((x,y)) 处的像素开始，把原始像素值替换为量化结果 $(\\bar{f}(x,y))$ ；阈值一般取 50%，高于阈值输出 1，低于阈值输出 0。随后计算量化误差：

$(\\delta=f(x,y)-\\bar{f}(x,y))$ 误差按如下规则扩散（见图 12‑19）：

将 $(\\displaystyle \\frac{3}{8}\\delta)$ 加到右邻像素；
将 $(\\displaystyle \\frac{3}{8}\\delta)$ 加到下邻像素；
将 $(\\displaystyle \\frac{1}{4}\\delta)$ 加到右下对角像素。

接着移动至同一行的下一个像素 ((x+1,y)) ，取更新后的值 $(f(x+1,y)+\\frac{3}{8}\\delta)$ ，重复上述量化‑误差扩散流程。一行处理完毕后再切换至下一行，循环往复。

边界处理规则：一行最右侧像素不再向右、向右下扩散误差；图像最后一行的像素不再向下、向右下扩散误差。

该抖动算法存在一个明显缺陷：误差会沿着对角线方向传递扩散，最终图像上会出现带有方向性的伪影。图 12‑20 展示了测试图像经过弗洛伊德‑斯坦伯格算法处理后的结果，可以清晰观察到对角线误差扩散带来的拖影效果。学术界已经提出了该算法的多种改进变体，用来消除这类方向伪影，更多内容可参阅 12.5 节。

彩色图像的抖动处理

前面主要介绍了灰度图像的 1‑bit 二值抖动；但同样的抖动技术也适用于彩色图像量化，可以避免或大幅减轻色带（量化轮廓）现象。例如图 12‑21 左侧图像直接从 24‑bit 量化压缩到 8‑bit，未开启抖动；右侧图像同样量化为 8‑bit，但是启用了弗洛伊德‑斯坦伯格抖动算法。对比可见，经过抖动处理的图像几乎看不到人眼可感知的色带痕迹。

12.4.2 基于空间填充曲线的抖动算法

正式展开讲解前，我们先回顾抖动算法的分类框架（图 12‑4），定位一下目前所介绍算法所处的类别，分类结果如图 12‑22 所示。

从分类图中可以看出：截至目前我们还没有介绍过非周期簇状抖动算法。理论上这类算法可以生成非常优秀的半色调效果，因为它能够模拟感光胶片乳剂的成像机理。在黑白胶片摄影当中，胶片表面涂满大量随机分布、互不重复的微小感光颗粒；每一颗感光颗粒都占据一块有限面积，曝光后整簇颗粒统一点亮或熄灭，最终生成照片。

图 12‑22 当中空缺的分类位置，就由一类重要的非周期数字半色调算法填补。这类算法使用一类分形曲线，统称为皮亚诺曲线（Peano curve），也常被称作空间填充曲线。相较于前面介绍过的算法，它具备多项优势：

通过计算每个抖动单元的平均量化误差并向外扩散，同时兼顾点簇聚合与网点离散分布两种特性；
平均量化误差扩散没有固定的传播方向，不会出现弗洛伊德‑斯坦伯格算法的方向性伪影；
输出属于非周期纹理，图像不会产生多余、规整的周期性摩尔纹路；
运行过程中可以动态修改抖动单元尺寸，让图像细节得到更好还原。

在数字出版领域，簇状聚合算法称为调幅（AM‑Amplitude Modulation）算法，原理是改变黑点簇的大小；网点离散扩散算法称为调频（FM‑Frequency Modulation）算法，原理是分散黑点并控制点的分布密度。本文介绍的空间填充曲线抖动算法效果介于调幅与调频半色调之间；因此可以灵活适配一大批网点定位精度各不相同的打印输出设备。

像素遍历排序

图像 $(f:U\\subset\\mathbb{R}^2\\to\\mathbb{C})$ 的像素遍历序列，是一条一对一映射路径 $(c:I\\subset\\mathbb{R}\\to U)$ 。该路径定义在实数域子集 I 上，路径扫过的区域 (c(I)) 覆盖图像定义域内所有像素。也就是说，遍历序列 c 会给坐标 ((i,j)) 的每一个像素分配唯一序号 k，满足 (c(k)=(i,j)) 。直观来讲，像素遍历就是按照一套固定顺序，不重复地走完图像当中所有像素。

矩阵存储格式图像最常用的遍历方式就是逐行扫描（见图 12‑23）。逐行遍历虽然简单高效，但用于抖动算法时有两个缺点：一是存在水平方向遍历倾向性；二是每行扫描结束跳转至下一行时路径出现断裂。往返式（牛耕式）遍历（见图 12‑24）虽然消除了换行路径断裂问题，但是遍历方向倾向性问题依旧存在。而正是这种方向倾向性，造成弗洛伊德‑斯坦伯格算法容易出现对角线伪影。我们可以采用近似空间填充曲线的遍历路径，一次性消除上述两个缺陷。

一条空间填充曲线是连续满射映射： $(\\phi:[0,1]\\to[0,1]\\times[0,1])$ 。该曲线可以连续遍历单位正方形内的每一个点。经过缩放变换后，就可以生成平面任意矩形区域上的空间填充曲线。我们的目标就是选取一条近似皮亚诺曲线的无重复路径，用来生成数字图像像素的遍历顺序。

空间填充曲线由数学家朱塞佩・皮亚诺于 19 世纪末首次发现。皮亚诺构造出了一条可以遍历单位正方形所有点的连续曲线，后人便将其命名为皮亚诺曲线；如今所有空间填充曲线也常泛称为皮亚诺曲线。图 12‑25 展示了希尔伯特曲线（Hilbert curve）递归构造过程当中的 4 个迭代步骤。不难看出：希尔伯特曲线每一次迭代生成的折线，都会恰好遍历规整网格上所有顶点，并且每个顶点仅访问一次。不断迭代即可生成任意分辨率图像对应的像素遍历序列；由此得到的遍历路径既不存在断点，也没有优先方向。图 12‑26 给出了 4 × 4 分辨率图像基于希尔伯特曲线生成的像素遍历顺序。

乍一看希尔伯特曲线遍历仅适用于分辨率为 $(2^p\\times2^p)$ （p 为正整数）的图像。不过任意 $(m\\times n)$ 尺寸图像对应的网格都可以嵌入一个边长为 2 的整数次幂的正方形网格；我们可以先对正方形网格生成希尔伯特遍历序列，再从中截取得到原图的遍历顺序。这种实现方案比直接构造矩形区域空间填充曲线遍历更加简便高效。除此之外还有许多其他递归构造的折线型空间填充曲线，均可用来生成离散图像像素遍历序列；下文示例统一选用希尔伯特曲线。

基于皮亚诺曲线的图像区域划分

想要理解皮亚诺曲线半色调算法，先要理解该遍历曲线如何将图像定义域分割成多个子区域。设 $(c:I\\to U\\subset\\mathbb{R}^2)$ 为基于皮亚诺曲线的图像像素遍历序列，将区间 I 拆分为 n 个子区间： $(I_1\\cup I_2\\cup\\cdots\\cup I_n)$ 。由于像素遍历属于一一映射，区间的分割会映射成图像区域的分割： $(c(I_1)\\cup c(I_2)\\cup\\cdots\\cup c(I_n))$ 。其中映射 $(c|_{I_j}:I_j\\to R_j=c(I_j))$ 即为子区域 (R_j) 内部像素的遍历序列。图 12‑27 展示一块 4 × 4 像素区块，被划分为像素数量分别为 5、4、7 的三块子区域（每个子区域使用不同灰度底色标注）。

值得强调：皮亚诺曲线生成的分割区块天然带有先后顺序。沿着皮亚诺曲线完整遍历一遍，就会按次序、不重复地访问每一块子区域。皮亚诺曲线划分出来的子区块在算法中的作用等价于有序抖动里的抖动单元，因此这些子区域同样被称为抖动单元，简称单元。

自适应分割

前文已经提到，非周期簇状抖动算法可以粗略模拟胶片乳剂半色调成像；而如果允许抖动单元大小动态变化，模拟效果还能进一步提升，这就对应胶片感光颗粒尺寸的变化。

图像中每一个像素都会归属某一个抖动单元，但一个抖动单元可包含多个像素，二者并不是一一对应关系。单元尺寸函数会为每一个像素分配一个数值，代表其所属抖动单元包含的像素总数。

想要在抖动过程中既保证灰度过渡自然，又精细还原图像细节，就应当依据像素邻域内亮度变化快慢动态调整抖动单元尺寸：像素周围亮度变化剧烈（高频细节区域）选用较小抖动单元；亮度平缓区域选用较大抖动单元。该变化量可以沿着空间填充曲线遍历方向，求取像素点亮度方向导数的绝对值得到。

实际上，单元尺寸与方向导数幅值之间的匹配关系，还要考虑到人眼对亮度变化的响应符合对数感知特性。因此单元尺寸函数取方向导数幅值的指数函数；该规则能够保证每个抖动单元内部的视觉亮度，与图像亮度沿着遍历曲线方向的变化斜率维持线性关系。图 12‑28 给出一幅灰度图像，以及采用希尔伯特曲线遍历生成的自适应分割效果图。

自适应分割计算流程

借助单元尺寸函数 s，就可以对图像定义域执行自适应区块划分，划分步骤如下：

设置用户给定的最大单元尺寸，作为当前单元尺寸初始值 (s_0)

；

按照空间填充曲线遍历顺序依次扫描图像内所有像素；

对当前像素 p，计算其目标单元尺寸 (s(p))

；

将

与当前单元尺寸

对比；若

，更新当前单元尺寸为 (s(p))

，然后继续扫描下一个像素 (p+1)

。

抖动单元量化

每个抖动单元的量化处理，需要依托单元内部像素的遍历顺序完成（见图 12‑29），详细流程：

使用式 (12.1) 计算抖动单元内部像素的平均亮度；

通过开启、关闭单元内部分像素，生成一张二值点阵图案，让该点阵的等效平均亮度尽可能接近单元原始平均亮度；二者之间的差值称为单元量化误差；

在生成二值点阵时，尽量将所有点亮的黑点聚集在一起，形成黑点簇。

图 12‑29 (a) 是一块包含 16 个像素、使用希尔伯特曲线遍历的抖动单元；图 12‑29 (b) 展示了该算法生成的黑点簇，分别对应该区域除去亮度 0 以外的 16 个灰度等级。图 12‑30 左侧是图 12‑28 当中灰度图像，每个抖动单元直接填充为区块平均灰度；右侧为采用上述量化方法生成的抖动效果图。

黑点簇位置微调

完成抖动单元量化后，我们会得到一块黑点簇。沿着空间填充曲线的走向滑动黑点簇，就拥有了一个位置调整的自由度。一种效果较好的策略：将黑点簇的中心像素对齐抖动单元的中心像素。该居中操作会为最终抖动图像引入一定随机性，从视觉观感上提升半色调画质，处理过程见图 12‑31：

(a) 给出一块 4 × 4 抖动单元；

(b) 定位单元内亮度最深的像素；

(d) 完成黑点簇居中放置。

完整抖动算法流程

简而言之，皮亚诺曲线抖动算法核心思路清晰。选定皮亚诺曲线完成像素遍历之后，依次执行下面步骤：

按照前文规则，将图像分割成一系列抖动单元；

使用式 (12.1)，计算每一块抖动单元内部图像的平均亮度；

为每一个抖动单元生成对应的黑点簇，保证生成二值图案的等效亮度近似等于单元平均亮度；

在抖动单元内部放置、居中对齐黑点簇；

沿着皮亚诺曲线前进至下一个抖动单元；计算下一块单元平均亮度时，叠加上一块单元遗留的量化误差；重复以上全部流程。

该算法执行效率很高，图像区块分割与各单元平均亮度计算可同步完成。图 12‑32 为测试图像经过皮亚诺曲线抖动处理后的结果；本例使用图 12‑25 中的希尔伯特曲线，抖动单元固定大小为 7 像素。图 12‑33 使用同一套算法，开启自适应单元模式，单元最大尺寸限定为 7 像素。

可以看到该算法沿着皮亚诺曲线的路径传递量化误差：当每个抖动单元仅包含单个像素，算法等价于沿着皮亚诺曲线做像素级误差扩散；当抖动单元尺寸更大时，则先在单元内部生成黑点簇，再将误差扩散传递给后续单元。

有序抖动当中，一旦选定抖动单元尺寸，点阵聚类图案就由抖动矩阵永久固定。而皮亚诺曲线抖动算法，黑点簇的尺寸完全取决于图像分割结果，程序运行期间就可以随时修改。

图 12‑34 截取了一幅卡通图像，画面包含线条轮廓以及大片灰度均匀区域；背景墙面阴影使用规则网点模拟传统印刷半色调网屏效果。图像分别经过固定单元尺寸、自适应单元尺寸两种皮亚诺抖动算法处理；自适应方案画质提升效果十分明显。改进的根源来自图像高频边缘区域的处理优化，自适应算法既能完美对齐线条轮廓边缘，又能利用均匀网点精准还原出各个灰度层次。

12.3 图案抖动